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EXECUTIVE  SUMMARY 


Energy  radiated  from  a  source  usually  propagates  in  many  directions.  To  detect  the 
spatial  origin  of  the  source  of  energy,  data  from  an  array  of  receiving  elements  (sensors)  are 
summed  coherently  (in  a  spatial  sense)  to  enhance  the  signal-to-noise  ratio  over  a  single 
element.  This  is  called  spatial  processing  or  beamforming.  Adaptive  beamforming  allows  for 
spatial  nulls  to  be  positioned  so  that  strong  spatial  interferers  can  be  cancelled.  In  a  multipath 
environment,  however,  multipath  signals  arrive  at  the  array  at  an  angle  spatially  different  from 
the  direct  path  as  well  as  being  time  delayed  (temporally  different)  from  the  direct  path. 

Unfortunately,  the  performance  of  most  adaptive  beamforming  algorithms 
progressively  deteriorates  with  the  amount  of  correlation  between  sources,  here  direct  and 
multipath.  In  addition,  most  adaptive  beamforming  algorithms  attempt  to  null  all  spatially 
different  sources  and  thus  cannot  make  use  of  the  multipath  information. 

This  report  introduces  methods  which  attempt  to  coherently  recombine  the  spatially 
and  temporally  different  multipaths  with  the  direct  path  to  further  enhance  the  array  gain. 
These  methods  are  based  on  utilizing  a  broadband  conventional  beamformer  incorporating 
adaptive  equalization  to  recombine  multipath  information.  Preliminary  results  are  compared 
with  two  well-known  narrowband,  purely  spatial  processing  techniques:  the  Bartlett 
beamformer  (also  called  a  frequency-domain  conventional  beamformer)  and  the  minimum 
variance  distortionless  response  beamformer.  While  preliminary  results  are  promising, 
several  important  research  issues  are  outlined  which  remain  to  be  addressed. 


gCtPX> 


111 


CONTENTS 

EXECUTIVE  SUMMARY  . iii 

SYMBOLS  . ix 

1.  INTRODUCTION . 1-1 

1.1.  Spatial  Processing  (Beamforming)  . 1-1 

1.2.  Temporal  Processing  (Channel  Equalization)  . 1-3 

1.3.  Signal  Model  . 1-3 

2.  GENERAL  BEAMFORMING  ISSUES  . 2-1 

2.1.  Performance  Versus  Array  Size . 2-1 

2.2.  Effects  of  Signal  Bandwidth  . 2-1 

2.3.  Source-to-Element  Propagation  Channel  . 2-10 

2.4.  Source-to-Element  Propagation  Channel  Model  . 2-12 

2.5.  Source-to-Element  Propagation  Channel  Equalization . 2-13 

3.  TIME-DOMAIN  BEAMFORMING  . 3-1 

3.1.  Continuous  Time-Domain  CBF  . 3-1 

3.2.  Discrete  Implementation  of  Time-Domain  CBF . 3-2 

4.  FREQUENCY-DOMAIN  BEAMFORMING  . 4-1 

4.1.  Frequency-Domain  CBF  (Bartlett) . 4-1 

4.2.  MVDR  . 4-2 

5.  BROADBAND  BEAMFORMING  INCORPORATING  ADAPTIVE 

EQUALIZATION . 5-1 

5.1.  Multichannel  Adaptive  Filtering . 5-2 

5.2.  Generalized  Sidelobe  Canceller  . 5-4 

5.3.  Broadband  CBF  Incorporating  Adaptive  Equalization  . 5-6 

6.  SIMULATIONS . 6-1 

6.1.  Three  Narrowband  Sources  with  Large  SNRs . 6-1 

6.2.  Single  Broadband  Source  with  Single  Multipath . 6-10 

7.  CONCLUSION  . 7-1 

8.  REFERENCES  . 8-1 

APPENDIX  A.  MATLAB  Routines 

APPENDIX  B.  Time-Domain  CBF  of  a  Linear  Array  with  Uniformly  Spaced  Elements 


v 


FIGURES 

1.1.  Relationship  between  some  popular  beamforming  algorithms . 1-2 

1.2.  Coordinate  system  . 1-4 

2.1.  MVDR  beampattem  for  single  CW  at  25.0Hz/- 14  deg;  FFT  sizes  of  32,  64, 256, 

and  1024  points . 2-3 

2.2.  MVDR  beampattem  for  single  CW  at  25.0  Hz/- 14  deg;  FFT  size  of  32  points  ...  2-5 

2.3.  MVDR  beampattem  for  single  CW  at  25.0  Hz/- 14  deg;  FFT  size  of 

1024  points  . 2-6 

2.4.  Array  pattern  for  5-element  ULA  pointed  at  - 14  degrees . 2-7 

2.5.  MVDR  beampattem  for  single  broadband  signal  incident  at  - 14  deg;  FFT  size  of 

1024  points  . 2-8 

2.6.  Adaptive  array  for  narrowband  signals  (data  and  weights  complex) . 2-9 

2.7.  Adaptive  array  for  broadband  signals  (data  and  weights  complex) . 2-9 

2.8.  Effects  of  a  single  discontinuity  in  the  propagation  medium . 2-10 

2.9.  Effects  of  a  single  discontinuity  in  the  propagation  medium . 2-11 

2.10  Channel  impulse  response  and  transfer  function  . 2-13 

2.11.  Inverse  channel  impulse  response  and  transfer  function . 2-14 

3.1.  Continuous  time  time-domain  conventional  beamforming . 3-1 

3.2.  Discrete  time  time-domain  conventional  beamforming . 3-3 

3.3.  Prebeamforming-interpolation  time-domain  beamforming . 3-4 

3.4.  Data  interpolation  of  a  single  element  for  time-domain  beamforming  . 3-4 

5.1.  Channel  models  for  two  sources  and  three  elements . 5-1 

5.2.  Multichannel  adaptive  filter  for  broadband  signals  . 5-3 

5.3.  Generalized  sidelobe  canceller . 5-4 

5.4.  Generalized  sidelobe  canceller  implementation  of  the  linear  constrained  adaptive 

beamformer  for  broadband  sources  [24]  . 5-5 

5.5.  Time-domain  CBF  with  adaptive  channel  equalization . 5-6 

6.1.  MATLAB  generation  file  for  first  simulation  . 6-3 


6.2.  TD-CBF  and  GSC  beamforming  for  three  narrowband  sources . 6-4 

6.3.  Bartlett  and  MVDR  beamforming  summed  from  5  to  45  Hz  for  three  narrowband 

sources . 6-4 

6.4.  TD-CBF  and  GSC  PSDs  along  - 14-,  0-,  and  + 30-degree  steering  directions  for 

three  narrowband  sources . 6-5 

6.5.  MVDR  beamforming  for  three  narrowband  sources;  FFT  size  of  1024  points  ....  6-6 

6.6.  MVDR  beamformer  output  for  three  narrowband  sources  along  fixed  steering 

directions  (top)  and  along  fixed  processing  frequencies  (bottom)  . 6-7 

6.7.  Bartlett  beamforming  for  three  narrowband  sources;  FFT  size  of  1024  points _ 6-8 

6.8.  Bartlett  beamformer  output  for  three  narrowband  sources  along  fixed  steering 

directions  (top)  and  along  fixed  processing  frequencies  (bottom)  . 6-9 

6.9.  MATLAB  generation  file  for  second  simulation . 6-12 

6.10.  PSD  of  transmitted  broadband  source  and  received  data  from  element  1  and  5  . .  6-13 

6.11.  TD-CBF  and  GSC  beamforming  for  a  single  multipath  . 6-14 

6.12.  Bartlett  and  MVDR  beamforming  summed  from  5  to  45  Hz  for  a  single 

multipath . 6-14 

6.13.  TD-CBF  and  GSC  PSDs  along  -14-  and  0-degree  steering  directions  for  a  single 

multipath . 6-15 

6.14.  MVDR  beamforming  for  a  single  multipath;  FFT  size  of  1024  points . 6-16 

6.15.  MVDR  beamformer  output  for  a  single  multipath  along  fixed  steering  directions 

(top)  and  along  fixed  processing  frequencies  (bottom)  . 6-17 

6.16.  Bartlett  beamforming  for  a  single  multipath;  FFT  size  of  1024  points . 6-18 

6.17.  Bartlett  beamformer  output  for  a  single  multipath  along  fixed  steering  directions 

(top)  and  along  fixed  processing  frequencies  (bottom)  . 6-19 

6.18.  Multichannel  adaptive  filter  weights  for  a  single  multipath:  presteering  at 

0  (top)  and  - 14  (bottom)  degrees  [Nj  =N2  =  75]  . 6-20 

6.19.  PSD  of  multichannel  adaptive  filter  output  for  a  single  multipath:  presteering 

of  0  degrees  [V/=A(2  =  75]  . 6-21 

6.20.  Weights  of  the  multiple  single-channel  adaptive  filter  for  a  single  multipath: 

presteering  of  0  (top)  and  -14  (bottom)  degrees  [Nj  =N2=75]  . 6-22 

6.21.  Weights  of  the  multiple  single-channel  adaptive  filter  for  a  single  multipath: 

presteering  of  0  degrees  [Nj  =150,  N2=0] . 6-23 

vii 


A.l.  Block  diagram  of  beamforming  signal  processing  computed  with  MATLAB  using 


various  sources  of  input  data .  A-2 

A. 2.  MATLAB  code  listings .  A-2 

B. l.  Single  source  ULVA .  B-l 

B.2.  ULVA  pattern  for  |  =  0.5,  <pm  =  0  degrees,  N  =  8,  fsample  =  =*  .  B-3 

B.3.  ULVA  pattern  for  j  =  0.5,  0  =  80  degrees,  M  =  8,  /samp/e  =  00  .  B-4 

B.4.  ULVA  pattern  for  j  =  0.5,  0  =  0  degrees,  M  =  8,  fsampie  =  250  Hz .  B-4 

B.5.  ULVA  pattern  for  j  =  0.5,  0  =  0  degrees,  M  =  8,  /Sflmp/e  =  1500  Hz  .  B-5 

B.6.  ULVA  pattern  for  j  =  0.25,  0  =  45  degrees,  M  =  8,  fsampie  =  00  .  B-5 

B.7.  ULVA  pattern  for  j  =  1.0,  0  =  45  degrees,  M  =  8,  fsample  -  00  .  B-6 

B.8.  ULVA  TD-CBF  averaged  beampattem  M  =  40,  d  =  30.0  m, 

0!  =  0  degrees,  =  25  //z,  SNR1  =  0  d®,  02  =  —  5  degrees, 

/2  =  45  Hz,  SNR2  =  -  5  dB,  fsample  =  750  Hz,  I  =  2 .  B-6 

TABLES 


2.1.  Direct  and  multipath  data  for  the  uniform  linear  vertical  array  in  figure  2.9 _ 2-12 


6.1.  Source  parameters  for  first  simulation . 6-1 

7.1.  Capabilities  of  beamforming  algorithms  . 7-1 


vm 


SYMBOLS 


A‘,Ar,At  . 

c  . 

d  . 

/  . 

M  . 

N  . 

pi  p 

ry,ry  . 

R(f)  =  E[X(f)X(f)H\  . 

R(u)  =  e[  X°(k,u)  X%k,u)H  ] 


U  . 

W . 

Xm(t )  . 

Xa(k,u)  ... 
y(k,u),  Y(f,u) 


magnitude  of  the  incident,  reflected,  and  transmitted  signals 

speed  of  propagating  wave  (m/s) 

array  element  spacing  (m) 

frequency  (Hz) 

number  of  elements 

number  of  FIR  filter  taps  or  weights 

instantaneous  array  output  power  and  average  array  output 
power 

spatial  spectral  correlation  matrix 

spatial-temporal  correlation  matrix 
(or  steered  covariance  matrix) 

number  of  time  samples 

required  steering  delay  a*  the  xm(t)  element  to  look  at 
direction  u 

array  steering  directional  unit  vector 

array  element  weight  vector 

received  signal  at  the  mth  array  element 

(includes  source,  interferes,  and  noise) 

aligned  (also  called  steered)  array  element  data  matrix 

time-domain  array  output  and  frequency-domain  output  for 
look  direction  u 
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1.  INTRODUCTION 


The  fundamental  aims  of  an  array  processing  system  [9]  are  to 

1.  detect  the  presence  of  propagating  energy 

2.  determine  the  location  of  the  source  of  propagating  energy 

3.  classify  the  source  from  its  radiation  waveform 

The  difficulty  in  designing  signal  processing  algorithms  to  achieve  these  fundamental  aims 
resides  in  the  nature  of  the  propagating  energy.  This  energy  is  composed  of  spatial  and  temporal 
components.  While  joint  spatial -temporal  processing  has  been  proposed  in  [15],  for  tract- 
ability,  the  assumption  is  usually  made  that  the  two  components  are  separable,  i.e., 
x(t,u)  =  y(t)  zfu).  Some  sequence  of  temporal  and  spatial  processing  algorithms  is  therefore 
developed  to  meet  the  design  objective.  A  typical  sequence  might  be  (1)  temporal  processing 
(element  data  conditioning,  i.e.,  automatic  gain  control,  low  pass  filter,  FFT),  (2)  temporal 
processing  (time  averaging  to  form  an  estimate  of  the  spatial  spectral  covariance  matrix),  and 
(3)  spatial  processing  (narrowband  beamforming).  Even  though  in  some  applications  the 
assumption  ot  spatial— temporal  separation  might  be  valid  (or  valid  enough),  there  are 
numerous  examples  when  such  an  assumption  is  not  valid.  For  instance,  when  an  attempt  is 
made  to  detect  a  quiet  source  in  a  multipath  environment,  useful  information  can  be  present  in 
the  multipaths.  However,  the  multipath  incident  energy  is  both  spatially  and  temporally 
different  from  the  direct  path.  In  such  a  case,  to  coherently  add  the  multipaths  to  the  direct 
path  signal,  spatial  processing  cannot  be  separated  from  temporal  processing. 

The  aim  of  this  report  is  (1)  to  clarify  the  distinction  between  spatial  and  temporal 
processing  in  array  processing  and  (2)  to  introduce  several  array  processing  algorithms  suitable 
for  real-time  applications  which  combine  both  spatial  anl  temporal  processing.  These 
algorithms  are  based  on  utilizing  a  broadband  conventional  time-domain  beamformer 
incorporating  adaptive  equalization.  This  is  a  preliminary  report.  The  merits  and  pitfalls  of 
these  algorithms  are  discussed  and  demonstrated  in  simulations  in  the  hope  of  generating 
ideas  and  improvements. 

1.1.  SPATIAL  PROCESSING  (BEAMFORMING) 

Beamforming  refers  to  spatial  signal  processing  algorithms  used  to  focus  an  array  of 
spatially  distributed  elements  (also  called  sensors)  to  increase  the  signal-to-interference- 
plus-noise  ratio  (SINR)  over  that  received  by  a  single  element  [1].  It  can  be  used  to  accomplish 
the  first  and  second  fundamental  aims  discussed  above.  Focusing  is  accomplished  by 
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coherently  (i.e.,  in-phase)  summing  the  element  outputs.  Coherency  is  achieved  by  correcting 
f  •'r  the  appropriate  phase  delay  in  each  element  (created  by  their  spatial  separation)  by  making 
reasonable  assumptions  about  the  propagating  wave’s  characteristics.  Any  a  priori  informa¬ 
tion,  such  as  sound  speed  profiles  or  source  range,  can  be  used  to  obtain  a  more  accurate 
estimate  of  the  phase  delay  between  elements  in  sensing  the  propagating  energy  wavefront. 

Beamforming  can  be  classified  in  two  broad  categories:  conventional  (or  nonadaptive) 
and  adaptive.  In  conventional  beamforming  (CBF),  each  element  output  of  the  array  is 
weighted,  delayed,  and  summed  to  align  an  incident  wavefront  coming  from  a  particular 
direction.  The  resulting  directivity  or  spatial  filtering  of  the  array  is  data  independent ;  that  is, 
the  weights  and  delays  are  predetermined  constants.  Adaptive  beanJorming,  on  the  other 
hand,  adjusts  its  weights  and  delays  to  the  observed  data.  For  instance,  a  null  can  be  steered 
automatically  in  the  direction  of  a  strong  interferer  to  suppress  its  effects  on  the  array  output.  If 
sufficient  observations  are  available,  adaptive  beamforming  can  usually  outperform  conven¬ 
tional  beamforming  in  achieving  its  fundamental  aims.  Only  linear  adaptive  algorithms  will  be 
discussed  in  this  text.  Nonlinear  adaptive  algorithms,  like  Volterra  filters  or  neural  networks, 
are  left  to  other  sources. 


Figure  1.1  shows  schematically  one  way  of  relating  several  of  the  more  popular  beam- 
for.uing  algorithms  (see  [31]  for  a  more  comprehensive  comparison  of  frequency-do!  nain 
narrowband  beamformers).  The  first  classifier  is  the  domain  in  which  the  beamforming 
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FREQ-DOMAIN 

(BLOCK) 


TIME-DOMAIN 

(ITERATIVE) 
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(BLOCK) 
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NARROWBAND  NARROWBAND  BROADBAND  NARROWBAND  BROADBAND 
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FROST  MULTI-  GENERAL-  FREQ- 
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CBF 

(DELAY  & 
SUM) 


Figure  1.1.  Relationship  between  some  popular  beamforming  algorithms. 


1-2 


algorithm  is  computed,  frequency  (via  the  DFT)  or  time.  The  block/iterative  classifier  refers  to 
whether  the  processing  is  computed  on  a  block  of  data  or  with  each  new  sample.  Narrowband 
beamforming  searches  for  the  direction  of  an  incident  signal  with  frequency  equal  to  that  at 
which  the  processing  is  being  performed.  Broadband  beamforming  searches  for  the  direction 
of  an  incident  energy  source  (sum  of  all  frequencies).  A  signal  is  usually  defined  as  narrowband 
when  the  bandwidth  of  the  incident  sources  is  much  less  than  the  reciprocal  of  the  propagation 
time  of  the  wavefront  across  the  array.  More  will  be  said  on  the  use  of  narrowband  beam¬ 
forming  with  broadband  data  in  Section  2.2. 

1.2.  TEMPORAL  PROCESSING  (CHANNEL  EQUALIZATION) 

Temporal  processing  refers  to  any  operation  which  can  be  performed  on  a  single 
channel  of  data,  say  from  a  single  element.  A  channel  includes  everything  between  the  trans¬ 
mitter,  i.e.,  the  source  of  energy,  and  the  receiver,  i.e.,  a  single  element  of  an  array.  The 
distinction  between  temporal  processing  and  spatial  processing  is  sometimes  confusing 
because  many  of  the  standard  temporal  processing  techniques  (anti-alias  filtering, 
quantization,  Fourier  transforms,  discrete  filtering,  etc.)  also  carry  over  into  spatial 
processing.  In  this  text,  we  will  be  primarily  concerned  with  the  temporal  processing  concept  of 
channel  equalization. 

Most  real-world  channels  are  nonideal  in  that  their  frequency  response  does  not  have 
constant  amplitude  and  linear  phase.  A  channel  equalizer  attempts  to  “invert”  the  frequency 
characteristics  of  the  communications  channel  so  that  the  frequency  response  of  the  cascaded 
channel-equalizer  combination  is  ideal  [12], [16],  When  the  channel  is  not  known  a  priori  or  is 
time  varying,  it  must  be  equalized  by  an  adaptive  equalizer.  An  adaptive  equalizer 
automatically  varies  its  impulse  response  in  accordance  with  the  temporal  variations  of  the 
channel  characteristics. 

In  a  multipath  channel,  the  received  signal  consists  of  the  sum  of  a  number  of  individual 
components  that  have  traversed  paths  of  differing  lengths  between  the  transmitter  and  the 
receiver.  Depending  on  the  phase  relationships  between  these  components,  they  can  either 
add  constructively  or  destructively  at  the  receiver.  The  channel  equalizer  must  sort  out  the 
multipaths  from  the  direct  path  and  coherently  recombine  all  components  into  a  single  output. 

13.  SIGNAL  MODEL 

In  this  text,  it  will  be  assumed  that  L  sources  (which  include  interferes  as  well  as  the 
target  source),  s^t)  l  =  0, 1,2,...,L  -  1,  each  generate  a  stochastic  plane-wav^  signal  which 
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propagates  in  an  isovelocity  medium  through  M  spatially  separated  array  elements, 
xm(t )  m  =  0, 1,2 -  1.  A  'ource  refers  to  any  incident  plane-wave  signal;  thus,  a  single 
discrete  energy  source  may  actually  generate  multiple  sources,  a  direct  path  plus  any  multi¬ 
paths  due  to  reflected  or  refracted  paths.  The  plane-wave  approximation  is  valid  in  an 
isovelocity  medium  for  sources  located  at  a  distance  from  the  array  equal  to  roughly  60  or  more 
times  the  aperture  of  the  array  (the  so  called  “far-field”  solution)  [9].  An  isovelocity  medium 
implies  that  a  wave  travels  at  the  same  speed  of  propagation  throughout  all  space.  Note  that  a 
far-field  source  can  generate  a  curved  wavefront  in  a  multivelocity  medium  (see  [40]  or  [18]).  It 
is  assumed  in  this  text  that  the  speed  of  propagation  and  the  exact  element  location  is  known. 

Each  source  sfjt)  produces  a  plane-wave  response  at  the  xm(t)  element  equal  to 


where 

xm  =  position  vector  for  the  xm(t)  element 

si  =  negative  directional  unit  vector  for  the  s£t)  source 

c  =  speed  of  propagation  of  the  acoustic  plane  wave 

nm(t)  =  noise  at  the  x m(t)  element  (assumed  to  be  spatially  uncorrelated) 

is  the  difference  between  the  arrival  time  of  the  plane  wave  from  s£t)  at  the  mth 

element,  xm(t),  and  the  arrival  time  of  the  same  wave  at  the  origin.  Thus,  for  xm  1  1,  or 
xm  =  (0, 0, 0),  this  time  difference  will  be  zero.  Figure  1.2  displays  the  coordinate  system  used 
throughout  this  text. 


(f)  =  ELEVATION 
ANGLE 

(+DOWN,  -  UP) 

0  =  AZIMUTH  ANGLE 
(OR  BEARING) 

(0  -  NORTH) 

x  =  r  cos0  cos 0 
y  =  r  cos 0  sin# 
z  =  r  sin0 


Figure  12.  Coordinate  system. 
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2.  GENERAL  BEAMFORMING  ISSUES 

This  section  discusses  several  adaptive  beamforming  issues  relevant  to  the  develop¬ 
ment  of  the  beamforming  structures  in  Section  5.  In  particular,  it  briefly  discusses  the 
performance  of  the  adaptive  array  with  respect  to  the  number  of  elements,  the  signal  band¬ 
width,  and  a  multipath  environment. 

2.1.PERFORMANCE  VERSUS  ARRAY  SIZE 

Array  gain  is  the  improvement  in  the  SINR  due  to  beamforming.  One  definition 
commonly  used  is  [9] 


array  gain  = 


SINRou 

SINR 


(in  Ch.  4.2.4  of  [9],  the  noise  term  implies  interferers  plus  noise  as  explicitly  stated  here;  see 
also  [8],[11]).  Note  that  if  a  large  interferer  is  cancelled  by  a  beamformer,  the  array  gain  given 
by  Eq.  (2.1)  can  be  quite  large.  The  array  gain  at  the  output  of  an  ideal  M  element  array  for  the 
case  of  no  interferers  and  spatially  white,  isotropic  noise  is  equal  to 

array  gain  =  10  log10(M)  dB  ^2  2) 

regardless  of  the  arrival  direction  (see,  for  example,  p.  37  of  [1 1]).  As  the  interference-to-noise 
ratio  (INR)  increases,  however,  the  array  gain  will  change  depending  on  (among  other  things) 
the  number  of  interferers,  the  direction  of  each  interferer,  and  the  power  in  each  interferer. 
For  instance,  if  a  10-element  array  effectively  cancels  a  30-dB  interferer,  then  the  array  gain 
given  by  Eq.  (2.1)  can  be  as  large  as  40  dB  [=  10  log(10)  dB  +  30  dB].  For  more  on  array  gain 
see  [5]. 

Any  adaptive  array  with  M  independently  controlled  elements  has  M  degrees  of 
freedom,  which  are  typically  distributed  as 

M  =  [no.  of  constraints]  +[no.  of  steer  directions]  +  [no.  of  interference  nulls]  (2.3) 
For  example,  the  minimum  variance  distortionless  response  (MVDR)  has  1  constraint  of 
constant  magnitude  in  the  single  steer  direction  and  so  has  a  capability  of  generating  Af-2 
interference  nulls  in  its  array  response.  When  more  interferers  are  incident  on  an  array  than 
the  array  is  capable  of  cancelling,  the  adaptive  weights  will  collapse  (to  zero)  as  the  incident 
INRs  increase  and  the  array  gain  will  deteriorate  [8]. 


22.  EFFECTS  OF  SOURCE  SIGNAL  BANDWIDTH 


Broadband  signals  create  (at  least)  three  separate  problems  in  array  processing.  The 
first  deals  with  the  fact  that  the  interelement  phase  shift  is  a  function  of  frequency.  To  see  this, 
consider  the  interelement  phase  shift  for  a  broadband  signal 

inter  -  element  phase  shift  —2  n  f  [rm(u)  -  rm_  j(u)] 
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where  t m(u)  =  x-  c—u  is  the  time  it  takes  a  plane-wave  coming  from  direction  u  to  propagate 
from  the  xm(t)  element  to  the  origin.  The  interelement  phase  shift  is  seen  to  be  linearly 
dependent  on  the  frequency  of  the  signal.  We  can  investigate  the  ramifications  of  this 
dependence  further  by  considering  the  uniform  linear  array  (ULA)  with  element  spacing  d. 
For  this  case  and  a  signal  incident  to  the  array  from  angle  0,  we  find 

inter  -  element  phase  shift  =  2  n  f  ^  sin(<p )  (ULA)  (2.4) 

From  Eq.  (2.4),  it  is  apparent  that  a  single  signal  withfsource  arriving  from  angle  <f>SOurce  can  be 
perceived  by  the  array  as  having  a  different  frequency  content,  farray,  and  arriving  from  a 
different  angle,  <parray,  as  long  as  the  interelement  phase  shift  is  unchanged,  or 
farray  sin(<parray)  -  f S0Urce  sin(4> source)-  This  can  happen  in  narrowband  beamforming  when¬ 
ever  fan-ay  *  f source  ■  In  these  cases,  the  source  appears  to  be  incident  from  an  angle 

<P array  =  sin-1  L  sin(<p source)  (ULA)  (2.5) 

Jsource 

For  narrowband  frequency-domain  beamforming,  two  situations  can  arise  where 
fanay  *  f source-  (1)  when  the  temporal  DFT performed  on  each  element’s  time  data  is  not  of 
adequate  size  so  that  spectral  leakage  between  frequency  bins  is  significant,  and  (2)  when  the 
spatial  spectral  correlation  matrix,  R(fj),  and  the  steering  vector,  E(f2,u)  (defined  in  Section 
4),  are  formed  at  different  frequencies. 

Figure  2.1  demonstrates  the  effect  of  an  inadequate  size  of  temporal  FFT,  with  the 
MVDR  adaptive  beamformer  processed  at  15.625  Hz  (bin  centered  for  a  1024-point  FFT). 
The  incident  signal  on  the  5-element  vertical  ULA  with  d=30-m  spacing  is  a  narrowband 
signal  at  25.0  Hz  coming  in  at  - 14  degrees  (c=  1500  m/s)  and  sampled  at  1000  Hz  (all  examples 
in  this  text  are  sampled  at  1000  Hz).  One  would  expect  the  MVDR  processor  to  yield  zero 
output  when  processing  is  done  at  any  frequency  other  than  25.0  Hz.  However,  figure  2.1  shows 
that  for  FFT  sizes  of  32, 64,  and  256  points,  the  spectral  leakage  from  the  25.0-Hz  source  into 
the  15.625-Hz  bin  creates  an  apparent  arrival  angle  of  -22.8  degrees,  as  described  by  Eq.  (25). 
As  the  frequency  bin  leakage  is  decreased  with  larger  FFTs,  the  desired  solution  is  obtained.  In 
figure  2.2,  the  MVDR  beam  pattern*  is  plotted  for  the  same  vertical  ULA  and  25.0  Hz  source 
for  fan-ay  -  1, 2, 3, ...,  60  Hz.  An  FFT  size  of  32  points  was  used  to  clearly  illuminate  the  arc  sin 
dependence  of  Eq.  (2.5)  (clearly  farray  f source  )■  Figure  2.3  plots  the  MVDR  beampattem 
results  for  an  FFT  size  of  1024  (with  50%  overlap).  The  figure  shows  the  correct  frequency  and 

*  The  beampattem  is  defined  in  this  report  as  the  array’s  output  energy  ( I  Ptf,  u)  df  )  as  the  direction  of 
look  is  varied.  It  is  a  real  valued  quantity.  J 
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Figure  2.1.  MVDR  beampattem  for  single  CW  at  25.0  Hz/- 14  deg; 
FFT  sizes  of  32,  64,  256,  and  1024  points. 


angle  of  the  incident  signal.  Unless  otherwise  stated,  a  1024-point  FFT  with  50%  overlap  will 
be  used  in  all  narrowband  frequency-domain  beamforming. 


The  second  problem  associated  with  broadband  signals  is  spatial  aliasing,  which  creates 
a  directional  ambiguity  in  the  beampattem.  Expressions  for  spatial  aliasing  are  based  on  the 
sampling  theorem  NO  TAG  and  are  a  function  of  the  incident  angle  of  the  signals.  The  worst 
case  for  spatial  aliasing  in  the  ULA  with  element  spacing*/  is  when  the  incident  angle  is  +/-  90 
degrees  (end-fire).  Under  these  conditions,  spatial  aliasing  will  occur  when 


d<  =  i 

c 


1 


2  fuAX  (Z6) 
where  /max  is  the  largest  temporal  frequency  of  the  incident  signals.  The  best  case  for  spatial 
aliasing  in  the  ULA  with  element  spacing  d  is  when  the  incident  angle  is  0  degrees  (broadside), 
where  spatial  aliasing  will  occur  whenever 


!<  = 


fmx 


(2.7) 


In  the  previous  example,  d  =  30  m,  c  =  1500  m/s;  we  can  expect  the  onset  of  spatial  aliasing  to 
appear  when  the  frequency  of  the  incident  signal  is  between  25.0  Hz  and  50.0  Hz,  depending  on 
the  incident  angle  of  the  signal  (see  Appendix  B  for  more  discussion  on  grating  lobes  and 
time-domain  conventional  beamforming  (TD-CBF)  with  a  ULA).  Spatial  aliasing  can  be 
viewed  clearly  as  grating  lobes  in  the  array  pattern.*  The  grating  lobes  can  be  clearly  seen  in 
figure  2.4,  where  the  array  pattern  of  the  5-element  vertical  ULA  is  plotted  for  a  steer  direction 


*  The  array  pattern  is  defined  in  this  report  as  the  complex  valued  weighting  given  to  an  incident  signal  as  its 
direction  is  varied  when  the  array’s  direction  of  look  is  fixed.  For  equally  spaced  elements,  it  reduces  to  the 
SDatial  Fourier  transform  of  the  arrav  weiehts. 
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of  - 14  degrees.  Figure  2.4  is  plotted  from  — 180  degrees  to  + 180  degrees  for  a  complete  view 
of  the  grating  lobe  structure.  Note  that  because  the  array  is  linear  along  the  Z  axis,  its  array 
pattern  is  azimuthally  independent.  Thus,  the  array  pattern  actually  forms  a  main  cone  and  a 
grating  cone  when  Eq.  (2.6)  is  satisfied.  It  is  straightforward  to  show  that  the  true  angle  (versus 
that  defined  by  the  grating  lobes)  for /  >  1/2 /max  lies  between 


Figure  2.5  plots  the  MVDR  beampattem  for  a  broadband  signal  incident  on  the  5-element 
vertical  ULA  from  — 14  degrees.  The  grating  lobes  are  readily  apparent  for  f  >  40  Hz. 

The  last  problem  associated  with  broadband  signals  is  that  the  array  pattern’s  main 
beam  width  is  a  function  of  frequency  (as  well  as  a  function  of  the  steer  direction).  This  can  also 
be  seen  in  figure  2.4,  where  the  main  beam  gets  wider  at  lower  frequencies,  and  thus  the 
resolvability  between  two  closely  spaced  sources  becomes  more  difficult  at  lower  frequencies. 
A  closed-form  approximation  for  the  main  beam  width  of  a  uniform  linear  array  is  found  in 
Appendix  B.  Broadband  beamformers  which  integrate  across  frequencies  (either  implicitly  or 
explicitly)  will  have  to  take  this  issue  into  effect,  especially  if  significant  amounts  of  energy  exist 
at  the  lower  frequencies.  The  simplest  method  is  to  filter  out  the  low  frequencies  (as  well  as  the 
high  frequencies  for  reducing  spatial  aliasing)  to  retain  the  required  beam  width  for  adequate 
resolution. 
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Figure  2.4.  Array  pattern  for  5-element  ULA  pointed  at  -14  degrees. 


2- 


JO 


ELEVATION  (deg) 


Figure  2.5.  MVDR  beampattem  for  single  broadband  signal  incident  at  -14  deg; 

FFT  size  of  1024  points. 


This  section  has  introduced  some  of  the  effects  of  signal  bandwidth  on  spatial  pro¬ 
cessing.  Insight  here  is  crucial  to  understanding  the  beampattem  output  and  the  operation  of 
various  beamformers.  Figures  2.6  and  2.7  show  simplified  adaptive  arrays  for  narrowband 
signals  and  broadband  signals,  respectively  [28].  A  single  complex  weight  is  capable  of 
supplying  the  required  gain  and  phase  adjustment  at  each  element  for  narrowband  signals,  as 
seen  in  figure  2.6.  The  tapped-delay  line  in  figure  2.7  permits  the  adjustment  of  the  gain  and 
phase  at  a  number  of  different  temporal  frequencies,  as  required  for  the  array  processing  of 
broadband  signals.  Signal  preconditioning  has  been  left  off  both  figures  for  clarity. 


Figure  2.6.  Adaptive  array  for  narrowband  signals  (data  and  weights  complex). 


Figure  2.7.  Adaptive  array  for  broadband  signals  (data  and  weights  complex). 
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23.  SOURCE-TO-ELEMENT  PROPAGATION  CHANNEL 


Two  identities  from  ray  theory  can  be  used  to  generate  a  more  complete  picture  of  the 
source  (including  interferences)  to  element  propagation  channel  in  a  nonhomogeneous 
medium.  Consider  figure  2.8  where  an  incident  wave  is  colliding  with  a  single  abrupt 
discontinuity.  We  find  that  [7] 

1.  Incident  angle  is  equal  to  the  negative  of  the  reflected  angle 

0,  =  ~<Pr  (2.9) 

2.  Snell’s  Law: 

c°s(<Pi)  _  cos(<pt) 

C  J  C2  (Z.1U) 

where  cx  and  c2  are  the  speed  of  propagation  of  a  plane  wave  in  medium  1  and  medium  2, 
respectively.  From  Eq.  (2.10)  it  can  be  seen  that  there  exists  a  critical  incident  angle  where  if 
0,  <  <pc  no  transmission  results.  The  critical  angle  is 


<pc  =  cos 


(2.11) 


£ 

The  ratio  of  the  sound  speeds  in  the  two  mediums,  is  called  the  index  of  refraction. 


REFRACTED  (TRANSMITTED 
WITH  CHANGE  IN  DIRECTION) 


Figure  2.8.  Effects  of  a  single  discontinuity  in  the  propagation  medium. 


The  magnitude  of  the  reflected  and  transmitted  waves  depends  on  the  characteristic 
acoustic  impedance,  Z  =  g  c ,  where  g  is  the  medium  density.  For  instance,  it  is  straight¬ 
forward  to  show  that  the  magnitude  of  the  reflected  wave  is  related  to  the  magnitude  of  the 
incident  wave  by  [7] 


=  Ai  \z2  sin(<Pi)  ~  z\  m*(0<)1 
|Z2  sin(<pi)  +  Zj  si/i(0,)J 

and  that  the  magnitude  of  the  transmitted  wave  is  [7] 

At  =  Ai  r  2  Z2  sinjtpj)  1 
[Z2  sin(<pi)  +  Zj  sin{<pt) J 


(2.12) 


(2.13) 
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Using  Eq.  (2.9) -(2. 13),  we  can  compute  the  magnitude  and  direction  of  the  reflected 
and  transmitted  waves.  For  example,  assume  that  a  source  is  2000  m  from  a  uniform  linear 
vertical  array  and  250  m  from  an  abrupt  impedance  boundary,  with  Z2  =  1.5/1.45  Z\  =  1.035 
Zi;  as  shown  in  figure  2.9.  Table  2.1  summarizes  the  multipath  and  direct  path  data  for  various 
element  depths.  The  data  show  quite  clearly  that  the  magnitude  of  the  multipath,  AT ,  is 
different  for  many  of  the  elements.  For  a  large  array,  each  element  will  most  likely  see  a 
different  propagation  channel.  In  fact,  it  is  possible  for  an  element  to  see  no  signal,  direct  or 
multipath,  from  the  source  or  to  see  only  a  transmitted  multipath  and  no  direct  path  (see,  for 
instance,  the  top  two  elements  in  Z2  in  figure  2.9).  Thus,  for  an  algorithm  to  use  the  multipath 
information  to  increase  the  array  gain,  each  signal-to-element  channel  must  be  treated 
individually. 

The  underlying  difficulty  with  including  the  multipath  data  is  that  it  is  temporally 
delayed  and  deformed  (due  to  nonideal  propagation  such  as  spreading,  absorption,  ducting 
and  nonideal  boundary  interactions)  as  well  as  spatially  different  from  the  direct  path.  Most 
adaptive  beamformers  treat  incident  energy  which  is  spatially  different  from  the  steer 
direction  as  interference  to  be  nulled  out.  This  is  equivalent  to  making  the  implicit  assumption 
of  all  source-to-element  propagation  channels  being  equal  and  no  multipaths  existing. 


Figure  2.9.  Effects  of  a  single  discontinuity  in  the  propagation  medium. 
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Table.  2.1.  Direct  and  multipath  data  for  the  uniform  linear  vertical  array  in  figure  2.9. 


Element 
Depth  (m) 

a 

4* 

(deg) 

4v  =  4* 
(deg) 

x  Reflect 
(m) 

Direct 
Path  (m) 

Multi- 
Path  (m) 

Diff 

(ms) 

190 

1.0 

— 

12.4 

863 

2000.9 

2047.8 

32.3 

220 

— 

13.2 

936 

2000.2 

2054.5 

37.4 

250 

1.0 

— 

14.1 

1000 

2000.0 

2061.6 

42.5 

280 

0.94 

0.5 

14.8 

1056 

2000.2 

2069.6 

48.0 

310 

0.52 

5.0 

15.6 

1107 

2000.9 

2076.9 

53.0 

340 

0.4 

WSEM 

16.4 

1153 

2002.0 

2085.2 

58.7 

370 

0.33 

8.9 

17.2 

1194 

2003.6 

2093.9 

64.8 

2.4.  SOURCE-TO-ELEMENT  PROPAGATION  CHANNEL  MODEL 


Each  source-to-element  propagation  channel  can  be  modeled  as  a  FIR  filter  with 
impulse  response 

i 

=  d(n  -  ntrect)  +  ^  Armi  d(n  -  r™‘“)  (2.14) 

i=  1 

for  the  mth  element  and  /  multipaths  [12].  By  aligning  the  received  element  data  to  point 
towards  the  direct  path,  we  have  a  channel  impulse  response 

/ 

hj.n)  =  6 <»)  +  A'mJ  6(n  -  +  n%m') 

1=1 

1 

=  6{n)  +  £  Armi  d(n  -  n°mi  )  (2.15) 

i  =  i 

We  can  take  the  Z  transform  of  Eq.  (2.15)  for  a  single  multipath  (/=1)  to  find  the  channel 
transfer  function  between  the  source  and  the  mth  element 

00 

Hm(z)  =  ^  hm{n)  z~n  =  1  +  .2  16) 

n  =  0  '  ' 

The  power  spectral  density  of  the  mth  element  is  now  related  to  the  power  spectral 
density  of  the  source  by 


PSDx.(f)  =  |//„C0|2  PSDAf)  (2.17) 

where 

\Hm(f)\2  =  1  +  (Arm)2  +  2  ATm  cos(2  n  f  i»®)  (2  lg) 

The  spectral  nulls  in  Eq.  (2.18)  are  characteristic  of  a  multipath  channel  and  can  be  seen  in 
figure  2.10  forAm  r-  0.8,  and  nm°  =  5.  It  is  seen  from  Eq.  (2.18)  that  the  depth  of  the  spectral 
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CHANNEL  IMPULSE  RESPONSE  CHANNEL  TRANSFER  FUNCTION 


0  5  10  15  20  0  fs/2 

TIME  DELAY  NORMALIZED  FREQUENCY 

Figure  2.10.  Channel  impulse  response  and  transfer  function. 

null  in  the  channel’s  transfer  function  increases  with  the  amplitude  of  the  multipath  (Am  r )  and 
that  the  distance  between  spectral  nulls  decreases  as  the  delay  (rim  °)  between  the  direct  and 
multipaths  increases. 

2.5.  SOURCE-TO-ELEMENT  PROPAGATION  CHANNEL  EQUALIZATION 

Even  though  the  multipath  can  be  distorted  by  the  propagation  channel  and  arrives  at 
the  array  from  a  direction  other  than  that  of  the  direct  path,  it  contains  energy  coherent  with 
the  signal  of  interest.  If  each  source-to-element  channel  can  be  equalized,  i.e.,  the  effects  of  the 
multipath  are  negated,  then  the  SINR  of  each  element  and  thus  the  array  can  be  enhanced.  (It 
should  be  noted  that  two  broadband  signals  are  fully  correlated  if  they  are  amplitude-scaled 
and  temporally  delayed  replicas,  while  two  narrowband  signals  are  fully  correlated  if  they  have 
a  fixed  phase  difference  relative  to  each  other.  Multipath  signals  are  usually  partially 
correlated  to  the  direct  path  due  to  nonideal  reflections  or  ’■^fractions.) 

Continuing  under  the  premise  that  the  inverse  filter  can  be  found,  we  can  find  that  the 
transfer  function  of  the  ideal  equalizer  or  inverse  filter  is 

(2.19) 

(2.20) 

(2.21) 


For  the  single  multipath  example  above 

1  +  2  j  =  0 

using  the  geometric  series.  Thus,  the  ideal  inverse  filter  impulse  response  is 

hmv(n)  -  Arm) 1  / for  i  =  0,  nrm,  2 
h%v(n)  =  0  otherwise 
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Note  that  hminv(n)  is  stable  if  and  only  if  \Amr  |  <  1.  Figure  2.11  plots  the  first  20  taps  of  the 
ideal  inverse  channel  impulse  response  (Eq.  (2.21))  and  the  ideal  channel  transfer  function.  It 
is  apparent  from  comparing  figures  2.10  and  2.11  that  the  cascaded  channel -inverse  filter 
results  in  a  distortionless  transmission  albeit  delayed  in  time.  This  is  the  aim  of  channel 
equalization. 

INVERSE  CHANNEL  INVERSE  CHANNEL 

IMPULSE  RESPONSE  TRANSFER  FUNCTION 


0  5  10  15  20  0  fs/2 

TIME  DELAY  NORMALIZED  FREQUENCY 


Figure  2.11.  Inverse  channel  impulse  response  and  transfer  function. 
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3.  TIME-DOMAIN  BEAMFORMING 

Time-domain  beamforming  techniques  are  used  in  communication,  in  transient  signal 
analysis,  or  for  wideband  sources.  This  section  briefly  introduces  time-domain  beamforming 
and  its  discrete  implementation.  It  introduces  notation  and  concepts  used  later  in  this  text. 

3.1.  CONTINUOUS  TIME-DOMAIN  CBF 

To  coherently  add  the  outputs  of  all  elements  in  the  time-domain  requires  the  proper 
delay  of  each  element  before  summing.  The  delay  is  determined  by  the  desired  “steer”  (or 
“look”)  direction.  The  summation  can  be  written  as 

M- 1 

y(t,u)  =  w’m  xm(t  -  t m(u)  )  (3.1) 

m  =0 

where 

u  =  steering  directional  unit  vector  for  the  ( 9,<p)th  direction 

— ♦  — * 

y  "1/ 

*m(u)  =  ~2Lc — =  required  steering  delay  at  the  xm(t)  element  to  look  at  direction  u 
wm  =  weight  (or  shading  coefficient)  for  the  xm(t)  element 
Figure  3.1  is  a  direct  implementation  of  Eq.  (3.1).  It  represents  pure  spatial  processing,  since 
the  temporal  delay  operation  has  a  distortionless  transfer  function  (constant  amplitude  and 
linear  phase). 


Figure  3.1.  Continuous  time  time-domain  conventional  beamforming. 

2 

The  instantaneous  power  in  the  beamformer  output,  P^t,  u)  =  \y(t,  u)|  ,  will  display  a 
peak  whenever  the  steer  direction  equals  a  source  direction.  (Recall  that 
ly(r,u)|2  =  y(t,u)  y'(t,u),  where  y\t,u)  is  the  complex  conjugate  of  y(t,u).)  In  low  SNR 
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environments,  Py{t,  u)  must  be  time-averaged  to  adequately  resolve  a  source  direction.  Since 
Eq.  (3.1)  is  a  summation  of  all  energy  coming  from  a  particular  steer  direction,  it  can  be  used  to 
detect  broadband  sources  as  well  as  narrowband  sources.  To  classify  the  source,  the  spectrum 
of  the  beamformer  output  must  be  analyzed. 

3.2.  DISCRETE  IMPLEMENTATION  OF  TIME-DOMAIN  CBF 

To  implement  any  algorithm  in  hardware  requires  that  the  continuous  input  data  be 
sampled  so  that  t  -*■  k  At,  where  k  is  a  positive  integer  and  At  =  —  is  the  sampling  interval. 

fs 

Sampling  Eq.  (3.1)  and  normalizing  each  iteration  to  the  sampling  interval  At,  we  find  the 
beamformer  output  can  be  written  as 

y(k,u)  =  WHX°(k,u)  (3.2) 

where  W  is  the  Mxl  weight  vector 

WT  =  [w0,  wj,  w2, ...,  wM_  x]  (3.3) 

and  Xa(k,u)  is  the  Mxl  perfectly  aligned  data  element  vector  (or  snapshot)  at  time  k 

Xa(k,u)T  =  [x0(k  -  r0(u)),  x^k  -  r j(u)),  x2(k  -  t 2(u)),  xM_1(k  -  rM_1(u))] 

=  [*£(*),  rf(k),  x°2{k),  ...,  xaM_1(k)]  (3.4) 

The  alignment  of  the  element  data  essentially  creates  a  broadside  steer  direction.  With  these 
vectors  defined,  the  instantaneous  power  in  the  beamformer  output  can  be  rewritten  as 

P^k,  u)  =  y(k,  u)  y'(k,  u)  =  WH  X°(k,  u)  Xa(k,  u)H  W  (3.5) 

Figure  3.2  shows  the  new  TD-CBF  structure  and  notation.  From  the  figure  it  can  be  seen  that 
time-domain  CBF  is  very  simple  to  implement.  Simply  retain  a  time-history  of  each  element’s 
input  data,  and  then  sum  the  properly  delayed  data  from  each  element  in  the  array.  The 
average  power  output  is  found  by  taking  the  expectation  of  both  sides  of  Eq.  (3.5)  to  give 

PJu)  =  WH  R{u)  W  (3.6) 

where 

R(u)  =  £[  X‘ (*,  u)  X‘(k,  5)"  ]  (3.7) 

is  the  spatial  temporal  correlation  matrix  [9]  (also  called  the  steered  covariance  matrix  in  [41]). 
If  the  input  process  is  ergodic,  then  a  time  average  will  approach  Eq.  (3.6)  (see,  for  example, 
Ch.  2  of  [14]). 

Unfortunately,  the  element  delays  are  rounded  to  the  nearest  integer  in  the  discrete 
implementation.  Thus,  the  discrete  time  element  delays  are  given  by 
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(3.8) 


rm(u)  —  round 


where  round[x]  rounds  x  to  the  nearest  integer.  (The  round  operation  is  inherent  in  the 
sampling  process.)  It  is  apparent  from  Eq.  (3.8)  that  the  steering  delays  will  be  rounded  to 
integer  multiples  of  the  sampling  period.  Thus,  the  higher  the  sampling  frequency,  the  more 
steer  directions  are  possible  and  the  better  the  delineation  of  the  array  beam  pattern. 


Often  interpolation  is  required  to  obtain  the  desired  delineation  of  the  array  pattern. 
Two  possible  techniques  for  data  interpolation  are  (1)  conventional  time  series  interpolation, 
which  effectively  resamples  each  element’s  data  at  a  higher  sampling  rate  [9]  and  (2)  data 
interpolation  by  a  linear  2-tap  FIR  filter  or  a  quadratic  3-tap  FIR  filter  to  interpolate  between 
data  values  while  retaining  the  same  sampling  rate  (see,  for  example,  Ch.  3.1  of  [4]). 

Figure  3.3  shows  the  most  intuitive  means  of  implementing  conventional  interpolation, 
called  low-pass  prebeamforming  interpolation  [26].  Each  element  is  resampled  at  a  higher 
sampling  frequency,  allowing  for  more  precise  control  over  the  steering  delays.  Enough  data 
memory  from  each  element  must  be  retained  to  account  for  +/-  maximum  steer  delay  (in 
number  of  samples),  which  is  a  function  of  the  resampled  sampling  frequency,  the  speed  of 
propagation,  the  element  location  (with  respect  to  the  origin),  and  the  steer  direction.  When 
the  number  of  steer  directions  is  greater  than  the  number  of  array  elements,  a  computationally 
more  efficient  method  of  interpolation  can  be  used,  called  low-pass  postbeamforming 
interpolation  [26].  The  interpolation  operation  can  be  separated  into  two  operations:  ( 1)1-1 
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zeros  added  between  each  successive  data  sample  and  (2)  interpolation  filtering.  /  determines 
the  amount  of  interpolation  and  is  assumed  to  be  an  integer  greater  than  1.  In  postbeam¬ 
forming  interpolation,  the  zero  adding  is  done  before  the  delay  and  sum  beamforming  com¬ 
putation,  while  the  interpolation  filtering  is  performed  after  the  delay  and  sum  beamforming. 
Both  techniques  result  in  similar  improvements  in  determining  the  angle  of  arrival  of  the 
source. 


Figure  3 3.  Prebeamforming- interpolation  time-domain  beamforming. 


A  second  method  of  data  interpolation  which  does  not  increase  the  data  rate  is  shown  in 
figure  3.4  for  a  single  element.  First,  the  rounded-down  integer  delay  is  used  for  coarse  data 
alignment.  Then  a  small  FIR  interpolation  filter  (either  2  taps  for  a  linear  interpolation  or  3 
taps  for  a  quadratic  interpolation)  is  used  to  compensate  for  the  difference  between  the  actual 
delay  and  the  rounded-down  integer  delay.  Essentially,  the  interpolation  filter  provides  for 
fine  steering.  In  most  cases,  the  generation  of  the  aligned  data  matrix  will  require  some  method 
of  data  interpolation  for  adequate  beampattem  delineation. 


DATA  ALIGNMENT  INTERPOLATION  FILTER 
(COARSE  STEERING)  (FINE  STEERING) 


Figure  3.4.  Data  interpolation  of  a  single  element  for  time-domain  beamforming. 


4.  FREQUENCY-DOMAIN  BEAMFORMING 


Frequency-domain  spatial  processing  is  an  especially  efficient  technique  when  the 
stochastic  signals  are  narrowband  and  sustained  rather  than  transient  [9].  This  section  briefly 
introduces  two  of  the  more  popular  frequency-domain  beamformers,  Bartlett  and  MVDR. 
References  are  provided  for  more  detailed  information. 

4.1.  FREQUENCY-DOMAIN  CBF  (BARTLETT) 

The  array  output  can  be  analyzed  in  the  frequency  domain  as  well  as  the  time  domain. 
In  the  frequency  domain,  time-delays  transform  into  frequency-dependent  phase  shifts.  Thus, 
taking  the  Fourier  Transform  of  Eq.  (3.1),  we  obtain 


M—  1 

(4.1) 

m  =  0 

where  Xm(f)  is  the  single  element  (also  called  along-channel)  temporal  Discrete  Fourier 
Transform  of  jtm(r)  in  a  discrete  implementation.  Eq.  (4.1)  can  be  rewritten  in  terms  of  linear 


matrixes 

r(f,Z) 

=  E(f,  u)H  W1*  X(f) 

(4.2) 

where 

E(f,u)T  =  [e+;'2*/r„©  , 

e+j2nfxl(u)  ,  *  e+;2n/T4f_,(5)j 

(4.3) 

WT  =  diag  [w0  ,  Wj  ,  .. 

•  »  wM-l\ 

(4.4) 

mT  =  [*o (f)  >  (/) , 

s 

H 

1 

(4.5) 

E(f,u)  is  the  plane-wave  steering  vector  which  corrects  for  the  phase  delays  at  frequency  / 
between  elements  when  the  array  is  steered  in  direction  u.  The  FFTs  are  usually  performed 
with  a  50% -75%  overlap,  weighted  with  a  good  FFT  window  (e.g.,  the  Hanning  window),  and 
are  of  an  adequate  length  (e.g.,  1024-8192)  for  good  frequency  resolution  NO  TAG.  Padding 
the  along-channel  data  with  zeros  for  good  frequency  resolution  is  equivalent  to  the  data  inter¬ 
polation  discussed  in  the  previous  section.  WT  is  a  diagonal  matrix  containing  the  element 
weights  (shading  coefficients).  For  unity  shading,  WT  reduces  to  the  identity  and  can  be 
dropped  out  of  the  equations.  The  instantaneous  power  in  the  array  output  is  given  by 

P&u)  =  Y(f,u)  r(f,u)  =  E(f,u)H  W"  X(f)  X(f)H  W  E(f,u)  (4.6) 

Or  taking  the  expected  value  of  both  sides,  we  find  the  average  power  output  as 
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where 


(4.7) 


PJf,u)  =  E(f,u)H  WH  R{f)  WE(f,u) 

R(f)  =  E[X(f)X(f)H] 

is  the  spatial  spectral  correlation  matrix  NO  TAG  (also  called  the  cross-spectral  density  matrix 
[6]).  In  practice,  only  an  estimate  of  R(f)  can  be  computed: 

N 

m-m-ji  (4.8) 

i  =  i 

The  efficiency  of  Eq.  (4.7)  is  that  once  an  estimate  of  R(f)  is  formed,  any  new  direction 
can  be  searched  by  simply  generating  a  new  plane-wave  steering  vector,  E(f).  For  broadband 
signals,  Eq.  (4.7)  is  computed  for  multiple  frequencies  within  the  band  of  interest  and  then 
summed  across  frequency. 


4.2.  MVDR 


MVDR  is  an  example  of  a  linearly  constrained  adaptive  beamforming  algorithm.  It  has 
the  ability  to  adapt  to  the  spatial  interferer  environment  and  place  spatial  nulls  wherever 
strong  interferes  exist.  The  MVDR  method  determines  the  weight  vector,  W(fu),  which 
minimizes  the  average  array  output  power: 

Ptf,u)  =  W(f,u)H  R(f)  W(f,u)  (4.9) 

subject  to  the  constraint  that  the  steer  direction  has  distortionless  response  (constant 
amplitude  and  linear  phase),  i.e.,  Wif,u)H  Eif,u)  =  1.  It  was  shown  in  [9]  (see  also  [31] -[40] 
for  more  on  MVDR)  that  the  array  element  weight  vector  satisfying  these  constraints  is 


wtfg>=  mum 

Ul  ’  E(f,u)H  R(f)-‘ E(f,u) 

Using  Eq.  (4.10)  in  Eq.  (4.9),  we  find  the  average  output  power  to  be 


(4.10) 


Eif, u)H  R(f) ~ 1  Eif, u) 


(4.11) 


Once  an  estimate  of  R(f)  is  formed,  its  inverse  must  be  computed.  This  requires  0(Af*) 
operations  in  general  and  can  be  a  numerically  unstable  operation  if  the  condition  number  of 
R(f)  is  large.  However,  once  this  inverse  is  computed,  any  direction  can  be  searched  by  just 
generating  a  new  plane-wave  steering  vector,  £($.  Since  the  spectra  of  MVDR  are  constrained 
in  the  direction  of  look  and  because  of  the  narrowband  processing,  MVDR  is  not  capable  of 
channel  equalization. 
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5.  BROADBAND  BEAMFORMING  INCORPORATING 
ADAPTIVE  EQUALIZATION 

This  section  discusses  adaptive  beamforming  algorithms  capable  of  channel 
equalization.  The  complexity  of  the  underlying  problem  is  depicted  in  figure  5.1  and  is 
described  in  more  depth  in  [13]  and  [15].  Each  source  signal  is  convolved  with  its  local  source 
impulse  response  (transducer,  ship  hull,  etc.),  Si(z),  before  being  transmitted  into  the  prop¬ 
agation  medium.  As  described  in  Sections  2.4.  and  2.5,  each  source-to-element  is  modeled  as  a 
separate  transmission  channel  (which  includes  multipaths,  distortion,  propagation  losses, 
etc.),  Hjj(z).  Upon  reception  at  each  array  element,  the  receiving  data  undergo  another 
convolution  with  the  element’s  impulse  response  (transducer),  Ej(z).  To  simplify  the  problem, 
it  will  be  assumed  that  each  element  is  equalized  and  that  the  problem  is  to  equalize  the 
channel  such  that  s,-  (t)  is  recovered.  The  channel  models  are  assumed  to  be  linear  and  slowly 
time  vaiying. 


CHANNEL  MODELS 


TRANSMITTING  PROPAGATION  MEDIUM  RECEIVING 

SOURCES  (MODELS  DEPEND  ON  SOURCE  ARRAY 

AND  ELEMENT  LOCATIONS) 

Figure  5.1.  Channel  models  for  two  sources  and  three  elements. 
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5.1.  MULTICHANNEL  ADAPTIVE  FILTERING 

Two  of  the  earliest  detailed  descriptions  of  the  multichannel  adaptive  filter  for  use  in 
adaptive  arrays  can  be  found  in  [28]  and  [23].  The  multichannel  adaptive  filter,  shown  in  figure 

5.2,  adapts  all  of  its  weights  to  minimize  the  mean  square  error  between  the  multichannel  filter 
output  and  a  desired  signal.  Figure  5.2  is  an  adaptive  filter  which  can  adapt  to  both  the 
temporal  and  spatial  characteristics  of  the  input  data  and  is  thus  capable  of  channel  equaliza¬ 
tion  and  interferer  cancellation.  The  weaknesses  of  this  multichannel  adaptive  filter  are  that 
(1)  source  localization  is  not  available  from  the  filtered  output  as  was  true  with  the  previous 
beamformers  (although  localization  information  might  by  using  the  converged  weights  values 
similar  to  the  techniques  used  by  frequency  ^omai.  inear  prediction  beamforming  [6]  and 
[9]),  and  (2)  it  requires  a  desired  signal,  dK  or  at  least  an  approximate  desired  signal. 
Although  these  characteristics  might  eliminate  the  multichannel  adaptive  filter  from  applica¬ 
tions  in  passive  surveillance,  a  desired  sequence  does  exist  in  communication  systems  (via  a 
training  sequence)  and  active  surveillance  (transmitted  pulse). 

To  investigate  the  multichannel  adaptive  filter,  consider  the  ith  (single)  channel  with  N 
=  N]  +  N2  +  1  taps.  The  time-domain  filtered  output  is  given  at  time  k  by 

yfk,u)  =  Wfkf  Xf(k)  (5.1) 

where 

W,<*)7  =  [w'i,_,y2(*)  ,  -  ,  wi0(k)  ,  ...  ,  H><VVi(*)]  (lxV)  weight  vector  (5.2) 

X°(k)T  =  x?(k  +  N2)  ,  ...  ,  *?(£)  ,  ...  ,  x?(k  -  Vj)]  (1  xN)  aligned  data  vector  (5.3) 

The  M-element  multichannel  filtered  output  is  the  sum  of  the  M  single-channel  filtered 


outputs,  or  in  matrix  formulation: 

m-  1 

y(k,u)  =  ^y,{k,u) 

1=0 

=  W(k)H  Xa(k)  (5.4) 

where 

=  [w0(*)7' ,  W,(k)T  ,  ... ,  j(fc)rj  (1  xMN)  vector  (5.5) 

Xa(k)T  =  [xa0(k)T  ,  X°(k)T  ,  ...  ,  XaM_l(k)T]  (1  xMN)  vector  (5.6) 

Written  in  this  manner,  it  is  easy  to  see  that  the  optimum  Wiener -Hopf  solution  which 
minimizes  the  MSE ,E[e(k)e*(k)J,  where 

e(k)  =  d(k)  —  y(k,  u)  (5.7) 

is  found  in  [28]  to  be 

Wopt  =  P(M)  (5.8) 
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or  in  [14]  for  detailed  single-channel  adaptive  filter  theory.  In  Eq.  (5.8),  R(u)  is  the  multi¬ 
channel,  spatial  temporal  correlation  matrix  and  P(u)  is  the  multichannel,  cross-correlation 
matrix  given  by 

R(u)  =  £[  Xa(k)  Xa{k)H  ]  (. MNxMN)  vector  (5.9) 

P(u)  =  E[  Xa(k)  d\k)  ]  (MNx  1)  vector  (5.10) 

Numerous  methods  of  solving  Eq.  (5.8)  without  actually  inverting  the  MNxMN  matrix 
(which  requires  0(M3iV3)  operations  in  general)  have  been  explored  in  the  past  30  years.  The 
recursive  least  squares  (RLS)  solution  (also  called  Kalman  Filter)  requires  0(M2N2) 
operations  [22].  While  the  shifting  property  does  not  exist  between  channels  as  would  be 
required  for  (super)  fast  RLS  solutions  (i.e.,  O(MN)),  it  does  exist  within  a  single  channel,  as 
shown  in  figure  5.2.  This  property  has  been  exploited  by  [19]  and  [25]  to  derive  two  different 
multichannel  RLS  lattice  algorithms  with  0(M2N).  The  application  of  stochastic  solutions  like 
least  mean  squares  (LMS)  to  the  problem  does  not  require  matrix  inversions.  This  results  in 
0(MN)  algorithms,  but  the  performance  of  these  algorithms  is  dependent  on  the  data  [28], 


The  data  alignment  helps  to  minimize  the  required  length  of  the  single-channel  FIR 
filter,  N.  By  interpreting  Eq.  (3.8),  we  see  that  the  length  of  the  filter  N  (as  well  as  other 
parameters)  determines  the  possible  steering  delays,  which  in  turn  determine  the  angles  at 
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which  the  multipath  can  arrive  and  still  be  equalized.  Thus,  by  aligning  the  data  toward  a 
particular  direction,  a  range  of  angles  is  defined  for  which  the  filter  can  equalize  multipaths. 
For  example,  the  uniform  linear  array  with  M  elements  can  equalize  multipaths  as  long  as  they 
arrive  within  -  sin-1 
steered  direction,  u. 

However,  since  the  adaptation  is  based  on  minimizing  the  output  error  power,  if  the 
same  desired  sequence  is  used  for  all  steer  directions,  the  weights  will  always  attempt  to 
equalize  the  same  signal  (as  long  as  it  is  in  the  range  of  possible  angles  defined  by  N  above).  The 
generalized  sidelobe  canceller  to  be  discussed  in  the  next  section  will  modify  figure  5.2  to  use  a 
TD-CBF  to  generate  the  desired  signal.  It  is  dependent  on  the  steer  direction  and  will  thus 
provide  more  directionality  information. 


f  3  fs\  *  <  +  sin't  2  d  T  i\  degrees  of  the 


5.2.  GENERALIZED  SIDELOBE  CANCELLER 

The  generalized  sidelobe  canceller  (GSC)  was  developed  in  [24]  where  it  was  shown 
that  the  linear  constrained  adaptive  beamformer  (like  the  Frost  algorithm  [21],  MVDR,  and 
[20]  and  [29])  can  be  implemented  as  a  combination  of  a  nonadaptive  beamformer  and  an 
unconstrained  adaptive  beamformer.  Figure  5.3  shows  a  block  diagram  of  the  GSC,  with  the 
spatial  filters  are  derived  from  an  array  of  elements.  The  GSC  is  essentially  a  spatial  version  of 
the  single-channel  adaptive  noise  canceller.  The  main  array  generates  the  desired  signal 
(determined  by  its  steer  direction),  which  can  be  corrupted  by  interferes  that  bleed  in  through 
the  sidelobes  of  the  nonadaptive  main  array.  Each  auxiliary  array  spatially  filters  one  interferer 
(determined  by  the  auxiliary  array’s  steer  direction);  then  its  output  is  adaptively  weighted  and 
subtracted  from  the  main  array  output,  resulting  in  the  signal  only  (ideally).  For  broadband 
interferes,  the  adaptive  weighting  must  consist  of  more  than  one  weight,  as  discussed  in 
Section  2.2. 
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For  the  GSC  to  keep  from  cancelling  the  signal,  it  is  imperative  that  the  auxiliary  arrays 
output  only  interferers;  that  is,  that  they  place  a  spatial-  spectral  null  in  the  direction  in  which 
the  main  array  is  steered.  For  fixed  arrays  with  regular  geometries  (like  the  ULA),  steer 
directions  can  be  determined  for  each  auxiliary  array  which  place  a  spatial  null  (of  the  auxiliary 
array’s  array  pattern)  exactly  in  the  direction  in  which  the  main  array  is  steered.  For  example, 
the  M— element  linear  auxiliary  array  will  have  a  spatial  null  in  the  <pmain  direction  when  the 
auxiliary  array  is  steered  in  a  direction  given  by  [27]: 

<Paux  =  sin~^sin(0mflJ  +  fori  =  integer  *  0  (5.11) 


Figure  5.4.  Generalized  sidelobe  canceller  implementation  of  the  linear  constrained 
adaptive  beamformer  for  broadband  sources  [24]. 


M- 1  different  directions  exist  which  satisfy  Eq.  (5.1 1).  In  addition,  [27]  presents  a  technique  for 
reusing  elements  of  the  main  array  to  form  the  auxiliary  arrays  with  minimal  GSC  degradation 
due  to  signal  correlation  between  the  auxiliary  and  main  arrays.  In  general,  however,  a 
data-independent  blocking  matrix  must  be  used,  as  seen  in  figure  5.4,  which  guarantees  that 
only  spatially  different  sources  (assumed  to  be  uncorrelated  interferers!)  are  passed  to  the 
multichannel  adaptive  filter.  Since  the  main  array  is  nonadaptive,  no  channel  equalization  is 
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possible  with  the  GSC.  The  performance  of  the  GSC  can  be  expected  to  deteriorate  in  a 
multipath  environment  because  the  multipath  allows  correlated  energy  to  enter  in  the 
auxiliary  array. 

53.  BROADBAND  CBF  INCORPORATING  ADAPTIVE  EQUALIZATION 

One  problem  with  the  multichannel  adaptive  filter  presented  in  Section  5.1.  is  that  for  a 
significant  number  of  elements,  the  equalization  tends  to  be  very  noisy  due  to  the  weight 
misadjustment  (a  function  of  the  total  number  of  weights,  MN).  By  treating  each  element  as  a 
separate,  single  channel,  as  in  figure  5.5,  the  number  of  weights  in  the  update  loop  can  be 
reduced  by  M  and  the  equalization  performance  enhanced  due  to  the  reduced  misadjustment 
noise,  but  now  the  spatial  cancellation  properties  have  been  reduced  to  that  of  a  nonadaptive 
TD-CBF.  This  technique  without  the  presteering  has  been  successfully  used  recently  in  an 
underwater  acoustic  communication  channel  [17].  In  communication  channels,  it  is  imperative 
to  equalize  each  element’s  data  before  summation.  Like  figure  5.2,  figure  5.5  relies  on  a 
desired  signal  which  could  be  supplied  by  a  TD-CBF  to  form  a  variation  of  the  GSC  presented 
in  the  previous  section. 


Figure  5.5.  Time-domain  CBF  with  adaptive  channel  equalization. 
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6.  SIMULATIONS 


Two  simulations  were  conducted  to  compare  the  performance  of  the  presented  beam¬ 
forming  algorithms.  The  first  simulation  considers  the  case  of  multiple  independent  narrow- 
band  sources  while  the  second  simulation  considers  a  single  broadband  source  in  a  single 
multipath  environment.  Both  simulations  consider  large  SNRs  only. 

6.1.  THREE  NARROWBAND  SOURCES  WITH  LARGE  SNRS 

The  first  simulation  compares  the  TD-CBF,  GSC,  Bartlett,  and  MVDR  beamformers 
for  the  case  of  three  separate  narrowband  sources  incident  on  a  10-element  linear  vertical 
array  with  large  SNRs.  The  MATLAB  generation  file  listed  in  figure  6.1  shows  the  details  of  the 
simulation  as  well  as  the  details  of  the  beamforming  and  plot  generation.  The  power  in  the 
three  narrowband  signals  differs  by  6  dB  and  all  sources  are  incident  at  different  angles,  as 
described  in  table  6.1.  The  array  is  “cut”  for  a  25-Hz  signal  (i.e.,  the  frequency  at  which  the 
distance  between  elements  is  half  the  incident  signal’s  wavelength).  Appendix  A  describes  the 
MATLAB  files  and  their  relationships  to  those  generated  in  [6].  The  additive  gaussian  noise 
has  a  variance  of  0.0001.  (The  adaptive  techniques  presented  here  will  actually  work  better 
with  more  random  noise  because  the  noise  helps  to  condition  the  data  covariance  matrix.) 


Table  6.1.  Source  parameters  for  first  simulation. 


Source 

Azimuth 

(deg) 

Elevation 

(deg) 

Frequency 
(Hz)  ’ 

Amplitude 

K3 

1 

90 

0.0 

25.0 

1 

0 

2 

90 

-14.0 

45.0 

2 

45 

3 

90 

30.0 

15.0 

4 

90 

Figure  6.2  summarizes  the  two  time-domain  beamforming  algorithms,  TD-CBF  and 
GSC,  while  figure  6.3  summarizes  the  two  frequency-domain  beamforming  algorithms, 
Bartlett  and  MVDR.  Figure  6.2  shows  the  effects  of  the  anti-(spatial)aliasing  filter  with  the 
TD-CBF.  Recall  that  any  incident  energy  at  a  frequency  greater  than  that  which  satisfies  the 
half-wavelength  spacing  of  the  array,  here  25  Hz,  shows  up  as  a  directional  ambiguity  due  to  the 
existence  of  grating  lobes.  When  not  completely  filtered  out,  the  45-Hz/-  14-degree  signal 
creates  havoc  in  the  output  beampattem.  The  antialiasing  filter  used  here  is  a  29-tap  band-pass 
filter  (5-40  Hz)  with  linear  phase.  The  apparent  weak  signal  incident  at  roughly  +60  degrees 
is  from  the  45-Hz/- 14-degree  signal  leaking  through  a  grating  lobe.  The  output  can  be 
improved  by  using  a  better  antialiasing  filter.  Figure  6.4  plots  the  power  spectral  density  (PSD) 
for  the  TD-CBF  and  GSC  output  when  they  are  pointed  in  the  three  directions  of  the  incident 
signals.  The  figure  clearly  shows  that,  as  expected,  the  GSC  is  better  able  to  spatially  filter  out 
the  signals  not  in  the  steering  direction. 


6-1 


Figures  6.5 -6.8  plot  the  beamforming  characteristics  for  the  MVDR  and  Bartlett 
frequency-domain  beamforming  algorithms.  A  grating  lobe  is  clearly  defined  by  the  45-H z/ 
—  14-degree  signal  in  figure  6.5.  The  15-Hz/+ 30-degree  signal  creates  the  sine  inverse  function 
in  figure  6.5  described  by  Eq.  (2.5).  It  is  due  to  spectral  leakage  (even  though  a  1024-point  FFT 
with  50%  overlap  is  used)  from  the  15-Hz/+ 30-degree  signal,  since  the  background  noise  is 
negligible.  Essentially,  MVDR  is  magnifying  this  minimal  energy.  When  the  same  simulation 
was  rerun  with  a  noise  variance  of  0.01,  MVDR  displayed  spectral  peaks,  without  the  sine 
inverse  effects,  as  desired. 

MVDR  seems  to  work  particularly  well  only  at  the  ideal  operating  frequency  of  the 
array.  For  instance,  even  though  the  15-Hz/+ 30-degree  signal  has  the  most  power  of  all 
incident  sources,  the  MVDR  algorithm  does  not  display  the  +30  degree  direction  with  the 
most  output  power  in  figure  6.6.  On  the  other  hand,  the  Bartlett  algorithm  gives  a  truer  output 
value  in  figure  6.8. 

Figure  6.3  plots  the  sum  of  frequencies  from  5-45  Hz  for  both  frequency-domain 
beamformers  for  comparison  with  the  time-domain  beamformer  outputs  plotted  in  figure  6.2. 
Note  the  similarities  between  the  TD-CBF  and  the  FD-CBF  outputs.  When  the  frequency- 
domain  beamformers  were  summed  over  5-25  Hz,  the  15-Hz/+ 30-degree  and  25-Hz/ 
0-degree  signals  are  displayed  cleanly  without  the  grating  lobe  interference  centered  at  +60 
degree  and  without  the  - 14  degree  peak. 
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%  NB_ULVA  **  *' 

%  Linear  array  of  10  elements  spaced  30  meters  apart  assumed  to  be 
%  arranged  along  the  z-axis  (vertical)  with  multiple  incident  narrowband 
%  sources 

%  Source  Azimuth  Frea  Elevation  Amp  Phase 


% 

1. 

90 

25 

0 

1.0 

0 

% 

2. 

90 

45 

-14 

2.0 

45 

% 

3. 

90 

15 

30 

4.0 

90 

%  Signal’s  frequency  (Hz) 
%  Signal’s  Amplitude 
%  Signal’s  Phase  (degrees) 

%  Signals  Azimuth  angle  (degrees) 

%  Signals  Elevation  angle  (degrees) 

%  Noise  vector 


%  number  of  signals 


%  conversion  from  degrees  to  radians 


%  initialize  parameters 
% - 

fa  =  1000;  %  Sample  increment  (1/1 000 Hz)  (sec) 

dt  *  l/(fa); 

T  -  4096;  %  Number  of  samples 

c  =  1500;  %  Speed  of  sound  (m/sec) 

M  =  10;  %  Number  of  elements 

%  Anay  of  element  positions  along  z  axis  (meters) 

P  =  (30)  •  (  zeros(l,M);zeros(l,M);0:(M-l)j; 

% - - — - 

%  generate  narrowband  complex  signals 
%  ■ 

f_sig  *  (  25  45  15]:  %  Signal’s  frequency  (Hz) 

amp_sig  ss  [  1.0  2.0  4.0];  %  Signal’s  Amplitude 

phase^sig  *  (  0.0  45.0  90.0];  %  Signal’s  Phase  (degrees) 

az_sijf  *  [  90  90  90];  %  Signals  Azimuth  angle  (degrees) 

eljsig  =  (  0  -14  30];  %  Signals  Elevation  angle  (degrees) 

v  *  0.01  *ones(  1  ,M);  %  Noise  vector 

%  Generate  complex  element  array  data 

dgr  =  (pi/ 180);  %  conversion  from  degrees  to  radii 

L  *  length(f  sig);  %  number  of  signals 

X=zeros(M,T); 
for  1=1  :L-\ 

u  «  [  cos(az_sig(l)*dgT)*cos(el_sig(l),dgr), ... 
sm(az_sig())*dgr)*cos(eI_sig(l)*dgr), ... 
sin(el_sig(l)*dgr)]; 

%  no  noise 

X  *  X+nb_sig_gen(f_sig(!),P.u^mp_sig(l),phase_sig(l).dVr.c,v*0); 
end  ~ 

u  i  [  cos(az_sig(L)*dgr)*cos(el_sig(L)*dgr),  ... 
sm(az_sig(L)*dgr)*cos(el_sig(L)*dgr), ... 
sin(el_sig(L)*dgr)]; 

%  including  noise 

X  =  X+nb  sig  gen(f  sig(L).P.u^mp  sig(L),phase  s!g(D,dt,T,c,v), 

%  - 

%  Beamforming 
% 

%  search  directions 

S&U 

%***•**•***•**•******•*••••••**••*•*•••••••••••**••*••••••••••*•***• 

%  freq-domain  beamforming 

. . . . . . . . 

%  beamforming  frequencies  must  be  exactly  equal  to 
%  temporal  FFT  discrete  frequencies 
F=(l:59)*fs/1024; 

Py«zeros(length(el).lenglh(F)); 
for  i*l:length(F) 

R  *  sfecsd(  X ,  F(i),  P  ,  c.  dt,  ’OverlapFFT:  1024’); 

PyB(:,i)=bfclass(  X,  R,  F(i),  az,  el.  P,  dt,  c); 

PyM(:,i)*invdr(  X,  R,  F(i),  az,  el,  P,  dt,  c); 
end 

%mesh(fl  ipudfdbCPy)’)) 

%contour(flipud(db(Py)’)) 

%pk>t(el,db((Py(:,26),Py(:,l5),Py(:,46)]))  %25.3, 14.6,44.9  Hz 
%plo«(F.db(  Py(77,:);Py(91,:);Py(121, :)])’)  %  -14deg.0deg30deg 
%ploJel,dl^(wm^^5:^))^*um(PyM(:,5:46,))jV)  %  sum  over  5-45Hz 

%  time-domain  beamformtng 

% . . . . . . 

N*  10; 

mu  *  0.0005;  %  normalized  1ms  parameter 

%  Linear  Phase  BPF  characteristics  required  for  TD  Beamforming 
%  (note  1 :  with  a  20  tap  FIR  filter,  the  0-5Hz  is  really  not  attenuated) 

%  (note  2:  group  delay  *  -  10  iterations ) 
a  bpf*[lj; 

btpf -  firl(29,[5/(fs/2)  4<V(fs/2)]);  %  fbegin  «  5  Hz  -  fstop  «  40  Hz 

for  m«l:M 

X(m,:)*fiUer(b_bpt*^bpf^C(m.:)); 
end 

Py  -zeros(len  gth(el),  1 ); 

B*toeplitz((  1  zeros(l,M-3)],{  1  -2  1  zeros(l.M-3)]), 
for  i»  1:30: 180 

%  conventional  time-domain  beamforming 
[PyC(i:if29),y]  -  id  cbf(Xju,el<i:i*29),P,dt,c); 

%  Generalized  Side k>bc  Canceller 
[PyG(i:i+29Xy]  *  gsc(X^z,el(i:i-f29),P,dt,c,N,mu,B); 
end 


Figure  6.1.  MATLAB  generation  file  for  first  simulation. 
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Figure  6.2.  TD-CBF  and  GSC  beamforming  for  three  narrowband  sources. 


ELEVATION  (deg) 

Figure  6.3.  Bartlett  and  MVDR  beamforming  summed  from  5  to  45  Hz  for  three 

narrowband  sources. 
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Figure  6.4.  TD-CBF  and  GSC  PSDs  along  -14-,  0-,  and  +30-degree  steering  directions 

for  three  narrowband  sources. 
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Figure  6.5.  MVDR  beamforming  for  three  narrowband  sources;  FFT  size  of  1024  points. 
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Figure  6.6.  MVDR  beamformer  output  for  three  narrowband  sources  along  fixed  steering 
directions  (top)  and  along  fixed  processing  frequencies  (bottom). 
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Figure  6.7.  Bartlett  beamforming  for  three  narrowband  sources;  FFT  size  of  1024  points. 
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Figure  6.8.  Bartlett  beamformer  output  for  three  narrowband  sources  along  fixed  steering 
directions  (top)  and  along  fixed  processing  frequencies  (bottom). 
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6.2.  SINGLE  BROADBAND  SOURCE  WITH  SINGLE  MULTIPATH 

The  second  simulation  compares  the  TD-CBF,  GSC,  Bartlett,  and  MVDR  beam- 
formers  and  the  multichannel  adaptive  filters  for  the  case  of  a  single  broadband  source 
incident  at  0  degrees  and  its  single  multipath  incident  at  -14  degrees  on  a  5-element  linear 
array.  The  MATLAB  generation  file  listed  in  figure  6.9  shows  the  details  of  the  simulation  as 
well  as  the  details  of  the  beamforming  and  plot  generation.  Figure  6.10  shows  the  PSD  of  the 
broadband  signal  and  the  received  data  at  element  1  and  element  5.  The  broadband  signal 
generated  by  the  MATLAB  file  “bb_sig_gen.m”  (see  Appendix  A)  is  composed  of  four 
separate  components:  (1)  15-25  Hz  band-pass  filtered  noise  which  is  amplitude-modulated  at 
2  Hz,  (2)  exponentially  low-pass  filtered  noise,  (3)  a  narrowband  signal  at  18  Hz,  and  (4)  three 
odd  harmonics  of  the  narrowband  signal  at  24.17  Hz.  The  effects  of  the  multipath  are  clearly 
seen  as  spectral  notches  (refer  to  Section  2.4)  in  the  received  data’s  PSDs.  The  channel  is 
modeled  after  the  simulation  outlined  in  Section  2.3.  The  FIR  filter  coefficients  for  each 
source-to-element  channel  are  defined  as 

b(l,:)=[  0.9  zeros(l,31)  0.9  0  0  0  0  0.0  0  0  0  0  0.0  0  0  0  0  0.00  0  0  0  0  0.00]; 
b(2,:)=[  0.9  zeros(l,3l)  0.0  0  0  0  0  0.9  0  0  0  0  0.0  0  0  0  0  0.00  0  0  0  0  0.00]; 
b(3,:)=[  0.9  zeros(l,31)  0.0  0  0  0  0  0.0  0  0  0  0  0.9  0  0  0  0  0.00  0  0  0  0  0.00]; 
b(4,:)=[  0.9  zeros(l,31)  0.0  0  0  0  0  0.0  0  0  0  0  0.0  0  0  0  0  0.85  0  0  0  0  0.00]; 
b(5,:)=[  0.9  zeros(l,31)  0.0  0  0  0  0  0.0  0  0  0  0  0.0  0  0  0  0  0.00  0  0  0  0  0.47]; 

which  defines  the  direct  path  incident  at  0  degrees  (interelement  delay  equals  0)  and  the  multi- 
path  incident  at  - 14  degrees  (interelement  delay  equals  5).  The  total  power  in  the  multipath  is 
roughly  0.8  times  that  in  the  direct  path. 

Figures  6.11  and  6.13  summarize  the  results  of  the  two  time-domain  beamformers,  and 
figures  6.12, 6.14-6.17  summarize  the  results  from  the  two  frequency-domain  beamformers. 
In  figure  6.11,  the  TD-CBF  displays  both  direct  and  multipath  incident  directions.  However, 
the  GSC  displays  no  multipath  at  - 14  degrees  because  the  direct  path  is  adaptively  filtered  by 
the  auxiliary  arrays  and  subtracted  from  the  main  array’s  multipath  output.  As  can  be  seen  in 
figure  6.13,  neither  beamformer  eliminates  the  spectral  notch  created  by  the  multipath. 

In  figures  6.14  and  6.15,  MVDR  is  seen  to  effectively  lock  in  on  the  24.17-Hz  narrow- 
band  signal,  giving  excellent  directional  data  at  24  Hz  for  both  the  direct  and  multipath  signals. 
However,  the  multipath  signal  is  displayed  at  a  power  level  significantly  lower  than  the  actual 
incident  power  level,  presumably  because  of  the  strong  correlation  with  the  direct  path  source. 
The  Bartlett  beamformer,  plotted  in  figures  6.16  and  6.17,  has  significant  difficulty  in  any 
directionality,  presumably  because  of  the  wide  beamwidths  of  the  5-element  array.  However, 
it  is  readily  apparent  from  the  PSDs  along  the  0-degree  and  -14-degree  direction  in  figure 
6.17,  that  the  two  signals  are  related.  This  was  not  the  case  for  the  MVDR  algorithm  in  figure 
6.15. 
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Figure  6.12  plots  the  su;  .ie  frequency-domain  beamformer  outputs  from  5-45  Hz 
for  comparison  with  the  time-domain  beamformer  outputs  in  figure  6.11.  The  MVDR  beam- 
former  displays  peaks  at  the  correct  directions  corresponding  to  the  incident  direct  and  multi- 
path  signals,  but  the  relative  power  between  the  signals  is  significantly  different  than  the  actual 
0.8.  It  is  interesting  how  much  better  the  TD-CBF  (using  the  same  antialiasing  filter  as  before) 
did  with  respect  to  the  Bartlett  beamformer.  One  might  have  expected  similar  results  between 
the  two,  as  was  the  case  for  narrowband  signals  in  the  previous  simulation.  This  might  be  due  to 
energy  at  frequencies  greater  than  45  Hz  or  to  differences  in  the  spatial-temporal  correlation 
matrix  and  the  spatial-spectral  correlation  matrix  as  suggested  by  [41]. 

The  weights  of  the  multichannel  adaptive  filter  described  in  Section  5.1  are  plotted  in 
figure  6.18  for  presteering  angles  of  0  degrees  and  — 14  degrees  and  N]  =A(?  =  75.  The  weights 
show  that  in  both  cases,  the  signal  with  more  power  (0-degree  incident  signal)  is  interpreted  as 
the  direct  path  by  the  filter.  Thus,  as  long  as  both  direct  path  and  multipath  lie  within  the 
multichannel  adaptive  filter’s  spatial  window  defined  by  N  (refer  to  Section  5.1),  the  output 
power  contains  no  directional  information.  However,  the  weights  could  be  further  processed 
to  obtain  directional  information.  Figure  6.19  shows  that  the  PSD  of  the  multichannel  adaptive 
filter  output  has  been  equalized  since  no  spatial  nulls  exist.  Figures  6.20  and  6.21  plot  the 
weights  of  the  multiple  single -channel  adaptive  filters  discussed  in  Section  5.3  for  Nj  =fy  =  75 
and  Nj  =150,  N2=0.  The  weights  are  seen  to  take  on  values  described  by  Eq.  2.20. 
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%  MPATH_ULVA 

%  This  script  file  generates  the  Uniform  Linear  Verical  Array 
%  element  time-series  inputs  for  a  multipath  environment.  The  multipath 
%  is  generated  by  propagating  a  signal  (broad  or  narrow  band)  through 
%  separate  channels  (modeled  as  an  ARMA  linear  filter)  between  each  sensor 
%  and  a  single  source. 

% 

%  written  by:  Rich  North  7-1-92 

. . . . . . 


%  initialize  parameters 

% - 

T  ■=  8192; 
fs  *  1000; 
dt  =  l/(fs); 
c  «  1450; 

M  a  5; 

X«=zeros(M,T); 

fO=24.17; 


%  Number  of  samples^*  13) 

%  Sample  frequency 
%  Sample  increment  (sec) 

%  Speed  of  sound  (m/sec) 

%  Number  of  sensors  in  array 
%  array  element  time-series  data 
%  primary  CW  wave  (d  *  lambda/2) 


%  Anay  of  sf  siggen  positions  along  z  axis  (meters) 

P  *  (30)  •  [  zeros(LM)^eros(  1  ,M);0:(M-1 )]; 

% - 

%  generate  signal 

% - - - 

%  (real)  broad  band  signal 
*»bbsig  gen(T,25-0. 1  xO.fO,  1 8,fs); 

% - 

%  propagate  signal  through  ARMA  channel 

% - 

rand( ’normal’); 
nvar  *  0.01; 
anoise  «  sqrt(nvar); 
u  *  zeros(l.T); 

%source-sen*or#l-M  channel 
Lfir*53; 

b»zeros(M,Liir); 

b(  1  ,:)=4  0.9  zcros(  1.3 1 )  0.9  0  0  0  0  0.0  0  0  0  0  0.0  0  0  0  0  0.0  0  0  0  0  01; 
b(2,  .)*=  0.9  zeros(  1,31 )  0.0  0  0  0  0  0.9  0  0  0  0  0.0  0  0  0  0  0.0  0  0  0  0  0  ; 
b(3,0«=  0.9zeros(131)  0.0000  00.0  0000  0.9  0000  0.0  0  00  00); 
b<4,.)=  0.9 zerosfUl)  0.0  0  0  0  0  0.0  0  0  0  0 0.0  0  0  0  0  0.85  00  0  0  0}; 
b(5,:)*=  0.9  zeros(  1.31)  0.0  0  0  0  0  0.0  0  0  0  0  0.0  0  0  0  0  0.0  0  0  0  0  0.47J; 
for  m*l:M 

u%filter(b(m,:),a4)+anoise*rand(u); 

X(m,:)*u; 

end 

% - 

%  Beamforming 

% - 

%  search  directions 

az-190]; 

el«(-90:89J; 

% . . . . 

%  freq-domain  beamforming 

%  beamforming  frequencies  must  be  exactly  equal  to 
%  temporal  FFT  discrete  frequencies 
%F*[  1 :59]*fs/l  024; 

%Py*zeros(length(eI),length(F)); 

%for  i*l:lengtb(F) 

%  R  *  sfecsd(  X .  F(»X  P ,  c,  dt,  ’OverlapFFT:  1024’); 

%  Py(:,i)«bfclass<  X,  R,  F(i),  az,  el,  P,  dt,  c); 

%  Py(:,i)«mvdt(  X,  R,  F(i),  az,  el,  P,  dt,  c); 

%end 

%plot(el,db((Py(  :,25),Py  ( :,  1 0),Py(  1 9))))  %24.4.9.8,18.5  Hz 
%pkH(F,db([Py(77..);Py(9l,;)J)’)  %  -14deg.0dcg 

%plot(eldb(sum(PyB(l:180,5:45)’)),el,db(sum(PyM(  1:180.5:45)'))) 

. . . . . . 

%  time-domain  beamforming 

%•••*•••••*••••••••••••••••••••••••••••••••••••••••••••••••••• 

%  LMS  filter  parameters 
N  ■  150; 

mu  *  0.05;  %  nJms  parameter 


mu  *  0.05;  '*>  nJms  parameter 

%  Linear  Phase  BPF  characterisncs  required  for  TD  Beamforming 


a  bpf  « (1); 

b  "bpf  «  fir  1(29, (5/(fs/2)  40/(fs/2)l);  %  ftegin  -  5  Hi 

for  m»l:M 

X(m,:)«fiJler(b_bpLa_bpf,X(m,:)); 

end 

9bPy*zeros(length(el),l); 

%B=toeplitz((  1  zeros(l,M-3))Jl  -2  1  zeros(l,M-3)]); 

%for  1*1:30:180 

%  conventional  time-domain  beamforming 
«(Py(»:i-f29),yJ « td_cbf(X,az.eI(i:i429),P,dw); 

%  conventional  time^domain  beamforming  w/  single  channel  equalization 
%(Py(i:i429Xy)  *  adeql(X4z^K»:i+29XP,dt,c3,mu4); 

%  Generalized  Side  lobe  Canceller 
%[Py(i:i+29),y]  -  gsc(X^eKi:i+29)J»dt,c.N.mu3); 

%  Multichannel  Adaptive  Filter 

%(Py(i:H29),y)  *  mch  af(X^z,el(i:i+29),P,di,c3inuA); 

%end 

%[Py(77),yl  -  gscfX^e!(T7U»dLc,Nanu3); 

%Pxxl  ■  spectrum^  y.  1024,75); 

%[*»y(91),y}  *  g*cPUz,el(91).P,dU,N,mu.B); 

%Pxx2  *  spectrum(y,  1024,75); 

%plot(F.db([Pxxl(l;59,l),;Pxx2(  1:59,1)')))  %  -14deg.0deg 


%  fbegin  *  5  Hz  -  fstop  *  40  Hz 


Figure  6.9.  MATLAB  generation  file  for  second  simulation. 
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Figure  6.11.  TD-CBF  and  GSC  beamforming  for  a  single  multipath. 
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Figure  6.12.  Bartlett  and  MVDR  beamforming  summed  from  5  to  45  Hz  for  a  single 

multipath. 
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Figure  6.13.  TD-CBF  and  GSC  PSDs  along  -14-  and  0-degree  steering  directions  for  a 

single  multipath. 
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Figure  6.14.  MVDR  beamforming  for  a  single  multipath;  FFT  size  of  1024  points. 
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FREQUENCY  (Hz)  ALONG  - 14  deg  AND  0  deg 

Figure  6.15.  MVDR  beamformer  output  for  a  single  multipath  along  fixed  steering 
directions  (top)  and  along  fixed  processing  frequencies  (bottom). 
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Figure  6.16.  Bartlett  beamforming  for  a  single  multipath;  FFT  size  of  1024  points. 
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Figure  6.18.  Multichannel  adaptive  filter  weights  for  a  single  multipath:  presteering  at  0 
(top)  and  - 14  (bottom)  degrees  [Nj  = 75]. 
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Figure  6.20.  Weights  of  the  multiple  single-channel  adaptive  filter  for  a  single  multipath: 
presteering  of  0  (top)  and  -14  (bottom)  degrees  [Nj  =  75J. 
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7.  CONCLUSION 


This  report  has  introduced  and  discussed  the  spatial -temporal  processing  capabilities 
of  several  adaptive  and  nonadaptive  beamformers.  Table  7.1  summarizes  the  results.  Only  the 
multichannel  adaptive  filter  was  found  to  be  able  to  adaptively  cancel  uncorrelated  interferes 
as  well  as  to  equalize  the  direct-path  source  with  its  correlated  multipaths,  true  spatial  - 
temporal  processing.  It  was  shown,  however,  that  some  a  priori  information  about  the  source 
signal  is  needed  for  equalization  and  that  the  multichannel  adaptive  filter  did  not  provide 
directionality  information  without  a  further  processing  of  the  weight  values.  The  frequency- 
domain  MVDR  and  the  time-domain  GSC  beamformers  were  both  shown  to  perform  only 
spatial  processing  (adaptive  interference  cancellation)  and  no  temporal  processing.  Further¬ 
more,  the  performance  of  both  algorithms  degrades  with  the  correlation  of  the  interferes, 
greatly  reducing  the  usefulness  of  these  algorithms  in  a  multipath  environment.  In  short, 
further  research  is  required  to  generate  a  true  spatial -temporal  beamformer.  One  possibility 
might  be  to  incorporate  blind  equalization  (equalization  without  a  desired  signal;  see 
NO  TAG)  into  the  multichannel  adaptive  filter  and  develop  a  method  of  extracting  directional 
information  from  the  adapted  weight  matrix. 


Table  7.1.  Capabilities  of  beamforming  algorithms. 


Algorithm 

Capabilities 

Adaptive  Interference 
Cancellation 

Source-Element 

Equalization 

Performance  Affected 
by  Source  Correlation 

Frequency  Domain 

Bartlett 

No 

No 

No 

MVDR 

Yes 

No 

Yes 

Time  Domain 

TD-CBF 

No 

No 

No 

Multi-Ch  AF 

Yes 

Yes 

No 

GSC 

Yes 

No 

Yes 

Multiple 

Single-Ch  AF 

No 

Yes 

No 
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Appendix  A 
MATLAB  ROUTINES 


This  Appendix  describes  the  MATLAB  routines  used  in  this  text  to  compare 
performance  between  various  beamforming  algorithms.  MATLAB  routines  are  patterned 
after  those  written  by  Jeff  Allen  [6j.  The  input-output  vectors  and  notation  have  been 
matched  to  those  used  by  Jeff  Allen  [6].  Figure  A.l  highlights  the  interaction  between  the 
MATLAB  programs.  Element  input  data  are  typically  stored  in  the  time-domain  arrays,  with 
T  time-domain  inputs  from  each  element  of  an  M-element  array  as 


T  time  samples 


■  x0(/i)  Xo(/i  -  1)  ... 

x^n  -  T  +  1)  ‘ 

J 

X  = 

x,(n)  x,(n  -  1)  ... 
x2(n)  x2(n  -  1)  ... 

x,(n  -  T  +  1)) 
x2{n  -  T  +  1) 

1 

xM-t (n)  xM-,(n  -  1)  ... 

xM.,(n  -T  +  1) 

M  elements 


The  Af-element  array  (X,YJZ)  locations  are  stored  in  the  array  P: 


(A.1) 


M  elements 


P  = 


x0  x,  x2 

y0  y  i  y2 

z0  2,  Z2 


y*i-i 

zM-\ 


(A.2) 


The  steering  matrix  can  be  formed  as  time  delays  for  multiple  azimuthal  angles,  0,  and  multiple 
elevation  angles,  (J>: 


n  az*n  el 


A  = 


T,(0o,<t>o)  Tt(0„</>o)  ...  ti(0o,0i)  7,(0,,^.,) 

r2(0(b0o)  T2(®1.0o)  —  T2(0O>0l)  r2(0J>0l) 

taOWo)  2^(0 u<p0)  ...  73(00,0,)  T 3(0„0,) 


-  T.(0(h0 2)  T,(0„02)  .. 

...  ~2[0(),<f>2)  T2(fil,<p2)  .. 

...  r3(0o,02)  r3(0,,02)  .., 


elements  (A-3) 

t 


or  phase  delays 


n  az*n  el 


A 


exp 


j  2  n  f 


>,(0o,0o)  fl(01.0o)  ...  *l(0<»0l)  T,(0„0,) 

22(.0(x<Po)  r2 (0„0o)  -  r2(0o ,0,)  T2(0„0,) 
*i(Qo<<Po)  2 2(0 \iipo)  —  7j(0o»0i)  r3(0,,0,) 


*i(0o.02)  r,(0„02) 
r 2(^0* 0 2)  r2(0,,02) 
r3(0o,02)  z $(0 \,<p 2) 


A  M  elements 

J  (A.4) 

t 


The  beamformer  output  power  is  stored  in  the  array  iy  (7  is  the  notation  used  in  [6])  for 
multiple  azimuthal  angles,  8,  and  multiple  elevation  ang’es,  4>, 


'MM  y{0 o,0i)  ><60,0 2)  K^o.03)  •••' 
,  _  ><^l.0o)  M^I.0.)  K®l.02)  X^I.03>  - 

r  y(^2.0o)  y(fi2*  0 1 )  y(d2,<t>2)  y(02,<Pi)  - 


(A.5) 


A-l 


Figure  A.l.  Block  diagram  of  beamforming  signal  processing  computed  with  MATLAB 

using  various  sources  of  input  data. 

The  following  MATLAB  programs  were  used  extensively  in  this  report  and  are  listed  in 
their  entirety  in  Figure  A.2. 


1.  steer.m  (modified  slightly  from  [6]) 

2.  bb  sig  gen.m 

3.  Xalign.c 

4.  td_cbf.m 

5.  mch_nlms.c 

6.  gsc.m 

7.  mch_af.m 

8.  arraypat.m 

9.  bfclass.ru  (modified  slightly  from  [6]) 

10.  mvdr.m  (modified  slightly  from  [6]) 

11.  steer.m  (modified  slightly  from  [6]) 


function  A  -  steer(  az .  el .  P .  c  .  ID .  string) 


%  STEER  -  creates  an  array  of  steenng  vectors  for  the  given  azimuth  and 
%  elevation  angles  as  given  in  [  1 J. 

% 

%  INPUT: 

%  a 2  -  vector  of  azimuth  angles  (deg) 

%  el  -  vector  of  elevations  angles  (deg) 

^  P  -  3xM  array  of  the  M  sensor  positions 
c  -  speed  of  signal  propagation 
%  fO  -  frequency  to  beam  form 

%  string  -  ’Phase*  output  phase  delay  (default  if  not  defined) 

%  Time*  output  time  delay  (seconds) 

% 

%  OUTPUT:  let  M  denote  the  number  of  sensor;  set 
%  n_az  ■  length(az) 

%  n~el  *  length(el) 

%  Then  the  output  array  A  is  M  x  (n_az  x  n_el )  where  each  column  has  the 
%  form 

%  T 

9r  string  -  ’Phase’  (default):  A(:.k)  •  (  exp(  i  2  pi  ID  P  u(k)  /  c  ) 

% 


%  string  *  Time’:  A(:.k)  ■  |  P  u(k)  I  c ) 

% 

%  where  u(k)  is  the  k-th  look  vector  of  length  M  determined  as  foUnwv- 
%  factor  k  as 
% 

%  k-p’n_»i  +  q  ( 0=  <  q  <  n_uz  ) 

9E 

%  Then  observing  the  MATLAB  convention  for  indices  starting  at  1.  we  get 
% 

%  |  cos(az(o+l))cos(el(p+l))  | 

%  u(k+1)«  |  sin(az(q+l))cos(el(p+l))  J 

%  j  sm(e((p+l))  | 

% 

%  REFERENCES: 

%\\\  Allen.  Jeffery  C. 

%  "Direction  Finding  Database  and  Techniques  (DFDT)** 

%  Code  635.  NOSC.  San  Diego  CA  92*52-5000.  December  1990 
% 

%  function  A  -  steer(  az ,  el ,  P .  c  ,  10 .  string) 


(Contd) 
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n_az  m  lcngth(az); 
n_el  *  length(el); 

u_a2  «  J  cos(az*pi/180);  sin(ar*pi/180);  ones(l.n  az)l; 
u^el  *  [  cos(el*pi/180);  cos(e]*pi/180);  sin(el*pi/T80)J; 

%  for  the  k-th  elevation  angle,  produce  the  look  vector  containing  all 
%  the  azimuth  angles 
A  *  zeros(3,  n_az*n_el), 
fork»l:n  el 

A(:,  (o_az#{k-l)+l):(n_az,k))  »  u_az  .*  (u_e!(:.k)  *  ones(  1  ,n_az)); 
end 


%  calculate  delay 
if  (nargin  <  6) 

A  «  exp(  i*2*pi'fO'P“A/c  ); 
else 

if  (strcmp(  string.  'Phase’ )) 

A  *  exp{  i,2‘pi'fO*P"A/c  ); 
elseif  strcmp{  string.  Time' ) 

A  =  P‘*A/c; 
else 

error(’ ...  siring  set  incorrectly  (Phase  or  Time)  ...*) 
end 
end 

bb  sig  gen.m 


function  s»bb^  si^jtenfT.nvarl  jivar2.f0.fljs) 

%  BB  SIG  GEN: 

%  This  matlab  function  generates  a  broad  band  signal  with  odd 
%  harmonics. 

% 

%  INPUTS: 

%  T  -  number  of  data  samples  (T>2*N1 
%  nvarl  -  variance  of  the  noise.  varfnHnlj 
%  nvar2  -  variance  of  the  noise.  var(n2(n)J 
%  ft)  -  primary  sinusoid  frequency  (+ harmonics) 

%  fl  -  secondary  sinusoid  frequency 
%  fs  -  sampling  frequency 

%  OUTPUTS: 

%  s  -  IxT  vector  of  broad  band  data 

% 

%  ALGORITHM: 

% 


% 

% 

% 

% 

°7c 

% 

% 

% 

% 

% 


I- 


Amp  Modulation 
c-fa(n) 


i 


n)(n) - 1  BPF  |-((n)-  • - +  — s(n) 


I _ I 


I 


I 


I . I  I 

o2(n)-~  |LPF|-l(n)--~ 


Sin(f0/fa)- 


*  I 

%  sin(3*ffVfs) - + 

*  I 

%  sin(5-f0/fs) - + 

%  sin(fl/fs) - - 1 

% 

% 


%  s-bb  sig  gcnfDtvarl.nvar2jPil.fe) 
%  wriite*n  by:  Richard  North  6/92 

. . . 


% - 

9e  generate  sinusoids 

% - 

for  i  *  1:2:5 

s  *  s  +  0.8  ^  i*cos(2*pi,i,f0/fs*t):  %  sinusoid 
end; 

phase  *  pi/4.6;  %  any  constant  phase 
$  «  s  +  u.6*cos(2*pi*fl/fs*i  +  phase);  %  sinusoid 

Xalign.c 


r . . . 

*  XALIGN: 

*  This  matlab  mex  function  delays  and  weights  the  ume  series 

*  data  from  multiple  sensors  and  stores  the  aligned  data  for 

*  more  processing. 

*  USAGE  FROM  MATLAB: 

*  Xa  =  Xalign(X.iau) 

*  TO  COMPILE: 

*  /matlab/bin/cmex  Xalign.c 

*  INPUTS: 

*  X  -  MxT  matrix  where  each  row  consists  of  the  time  series 

*  data  from  a  single  element 

*  M  -  number  of  elements 

*  T  -  number  of  data  sampels 

*  TAU  -  Mxl  matrix  of  delays 

*  OUTPUTS: 

*  Xa  -  MxT  matrix  of  weighted  and  delayed  X  elements 

*  ALGORITHM: 

*  «  X(m,t-tau(m,u)) 

*  Xa(m.t) 

*  =  0  (if  X(m,t~  tau(m.u))  not  in  matrix) 

*  REFERENCE: 

*  R  North,  "Broadband  Conventional  Beamforming  Incorporating 

*  Adaptive  Equalization,"  NRaD  Tech.  Note  1992 

*  written  by:  Richard  North  5/92 

. * . . . * . •••••/ 

#include  <math.b> 

♦include  <stdio.b> 

♦include  "craexJT 

r  input  arguments  7 
♦define  X  IN  prhsfO] 

♦define  TAU  IN  prhs[l] 

r  output  arguments  7 
♦define  Xafigned_OUT  plhs(0] 

user  fcn(nlhs,plhsjirhs,prhs) 
int  nlhsjirhs; 

Matnx  •plhsfl,  *prhs(J; 

/•  Matlab  variables  7 
unsigned  int  M.Nlooks.T; 

double  *  X_re,*X_nn,*TAU,*  Xaligned^rc,*  Xaligncd_im; 

/•  misc  local  variables  7 
int  un.n^dc  lay  .delay _a; 

r  check  number  of  arguments  7 

if  ( nrhs !  *  2)  mex_error( " Xalign  requires  2  inputs."!; 

if  (nlhs  f ■  I )  mexJerTor("Xaiign  requires  I  output."); 


% - 

%  initialize  variables 

% - 

t  «  0:1  :T—  1 ;  %  t  -  0.1.2J.-..T-1  time 
anouc)  -  sqrtt nvarl); 
a  none  2  -  sqrt(nvar2); 
s  *  zerosf  1.1}; 
rand( 'normal ); 


r  set  variables  and  allocate  memory  7 
M  ■  TAU  IN-  >m;  /•  ♦  of  rows  equals  ♦  of  elements  7 

if  (X_IN->m  !-  M) 

mex_error( "Xal ign  requires  X,  and  tau  to  have  same  number  of  elements."); 
Nloolu  ■  TAU  IN- >n;  f*  ♦  erf  columns  equals  ♦  of  steer  directions  7 
if  (Nlooks  I- IT 

mex  error("Xalign  requires  only  1  look  direction."); 

T  *  )C  IN- >n;  /*  ♦  of  columns  equals  ♦  of  time  samples  7 


% - 

%  generate  amplitude  modulated  bandpass  noise 


n  -  &noisel*rand(i);  %  random  noise 
flow  -  15/(fs/2); 
fhigh  -  25/(fs/2); 

fb^ j»butteT(5/flow  fhighl);  %  1 0th  ordered  band  pass  filter 
m  *  filter(b  jui);  %  band  pass  filtered  nouc 
C  -  0.5;  %  C«0.0  -  >  no  mod 

s  *  (oncs(t)  +C*cos(2"pi*2/fs*t))."fn;  %  amp  mod  (2  H2)  bpf  noise 


% - 

%  generate  exponentially  filtered  low  pass  noise 


n  *  anoise2*rand(t); 
f  -  1 0:10 1/10; 


b  random  noise 


. -(0:10V5)  J; 

b*fir"2(  ibijn);  %  10th  ordered  low  pass  filter 

•-(1); 

t  ■  s  +  filter(bAn);  %  signal  +  filtered  noise 


X  re  -  X  IN->pr, 

X  im  »  >T  IN-  >pi; 
TAU  -  TXUJN-  >pr. 


if  (X_im  --  NULL)  r  real  7 

Xaligned  OUT  -  create_macrix(M,TREAL). 
Xaligned_re  ■  Xahgned_OUT->pr, 

else  r  complex  7 

Xahgned^OUT  ■  create  matrix(M.TCOMPLEX); 
Xaligned_re  *  Xaligned  "OUT -  >pr, 

Xaligned“im  *  Xaiignetf  OUT-  >pi; 


if  (X  im  --  NULL)  r  real  7 
< 

/*  do  delay  and  weighting  for  each  time  t  */ 

r  NOTE:  Matlab  stores  by  columns  so  that  X(u)*X[i+TbtalNrows*j] 
If  we  only  have  vectors,  this  » transparent  7 


Figure  AJ..  Continued. 
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for  (t*0;t<T;t+  +  ) 

/•  across  elements  •/ 
for  (m*0;  m<M;  m+  +  ) 

i 

a  m  t-(mt)TAU[m]; 

if  ((a  >=  0)  &&  (a<T)) 
Xaiigncd_re|m+M*t]  = 
else 

Xahgned  re(m+M*tj  «= 
}  /*  end  m  ’/ 


X_re[m+M*aj; 

0; 


}  /*  end  i  V 

) 

else  /•  complex  •/ 

/•  do  delay  and  weighting  for  each  time  t  */ 

/•  NOTH:  Matlab  stores  by  columns  so  that  X(ij)=Xfi+TbtalNrows*jj 
If  we  onlv  have  vectors,  this  is  transparent  •/ 
for  (t*0;  t<T;  t+  +  ) 

< 

/*  across  elements  •/ 
for  (m*0;  m<M,  m  +  +  ) 

a  *  t-(int)TAU(m]; 

delay  a  *  m+M*t;  i*  new  matrix's  location  •/ 
if  f  (a">  *  0)  &&  (a  <  T)) 

delay  *  m+M’a;  /•  location  within  old  matrix  • 
XaJ(gned_re|delay_a)  *  X  re[delayj; 
Xaligned3im(delay]_a]  =  Xjm|  delay); 

else 

{ 

Xaligned_rc(delay_a)  *  0; 

Xaligned_im(delay~  a)  *  0; 

}  /•  end  m  V 

}  f*  end  l  V 

} 


td_cbf.m 


function  [Py.Y]  =  td_cbf(  X,  az,  el.  P.  dt.  c) 

%  TD  CBF  -  classical  time  domain  "delay  and  sum"  beamforming 

% 


%  INPUT: 

%  X  -  MxT  array  of  T  time  samples  of  the  M  sensor  outputs 
%  a z  -  beamforming  azimuth  angles  (degrees) 

%  el  -  beamformine  elevation  angles  (degrees) 

%  P  -  3xM  array  of  the  M  sensor  positions  (ft)  (m) 

9c  dt  -  sample  increment  (sec) 

%  c  -  speed  of  sound  in  propagation  media  (ft/sec)  (m/sec) 


%  OUTPUT: 

%  Py  -  array  of  classical  beamformer  power  output;  Py(k,l)  is  the  output 
%  at  azimuth  angle  az(l)  and  elevation  angle  el(k). 

%  Y  -  NlooksxT  array  of  T  time  samples  of  the  beamformer  output 


%  ALGORITHM: 

%  As  discussed  in  [1],  the  time  domain  classical  beamformer  has  output 
%  H 

%  y(t)  *  [sum  (m*l)  ~(M>  X(m:t-tau_m))  (sum  (m«l)  ~  (M>  X(m:t-tau  m)J 

% 

%  for  steering  delay  vector 
%  H 

%  a  *  [p_(m>  u/c))  (seconds) 

% 

%  where  p_(m)  is  the  position  vector  of  the  mth  sensor,  u  is  the  look 
%  direction. 

% 

%  REMARKS: 

%  (1)  to  beam  form  (ie  delay)  in  a  arbitrary  desired  direction  requires 
%  interpolation  due  to  the  discrete  nature  of  the  input  time  samples. 

%  (  For  other  techniques  see  R.  Pridham  and  R.  Mucci,  Digital 

%  Interpolation  Beamforming  for  Low-Pass  and  Bandpass  Signals " 

%  Proc.  of  IEEE.  vol.  67,  no.  6,  pp.  904-919,  June  1979. ) 

%  (2)  y’  is  the  conjugate  transpose  of  the  vector  y 
% 

%  REFERENCES: 

%  [11 D,H.  Johnson  and  D.E.  Dudgeon. 

%  Fundamentals  of  Array  Signal  Processing", 

%  Prentice- Hall  to  be  published  in  1993. 


Y  ■  zeros!  N5ooks.T); 
Py  «  zeros(l.Nlooks); 
x  =  zeros(Tl ); 


9  Generate  steering  time  delays  for  ail  directions  of  look 

A  =  steer(  az .  el,  P ,  c ,  1 ,  Time');  %  delays  in  seconds  (freq  not  used) 
TAU  *  round(A“fs);  9c  delays  in  iterations 


9c  align  data  matrix  and  equalize  channels  for  each  look  direction 
for  n»  I  :N looks 

9c  align  all  clement  data  for  look  direction  n 
Xa  =  xalign(X.TAU(:ji)); 

9c  weight  all  element  data  for  look  direction  n 
9cX a  *  diag(W)*Xa; 

9c  equalize  each  element's  channel  i 
Y(n.:)  *  sum(Xa); 
end; 

Py  =  mean{Y’.*conj(Y’)), 

9c  tidy  up  and  reshape  the  output; 

9c  we  know  Py  must  be  real ... 

Py  *  rcal(Py); 

%  ...  and  return  y  with  the  rows  containing  the  azimuth  for  a  fixed 
%  elevation  angfe. 

Py  *  reshape!  Py  •  n_az .  n_el )’; 

mchjnims.c 


. . * . 

•  MCHNLMS.C: 

•  This  matlab  mex  function  uses  the  normalized  LMS  algorithm  to 

•  filter  the  multiple  channel  reference  input  X  to  generate  the 

•  filtered  output  y.  y  is  subtracted  from  the  primary  input,  d 

•  to  generate  the  error  output,  err. 

•  USAGE  FROM  MATLAB: 

•  [y,erT,wts}= mch_n]ms(X,d.N  jnu^lpha) 

•  TO  COMPILE: 

•  /matlab/bin/cmex  mch_nlms.c 

•  INPUTS: 

•  X  -  MxT  vector  of  reference  inputs 

•  T  -  number  of  data  samples 

•  d  -  lxT  vector  of  primary  input 

•  N  -  number  of  filter  weights/per  channel 

•  mu  -  LMS  step  size 

•  alpha  -  leak  parameter 

•  (alpha* 0  results  in  "standard  LMS" 

•  l/(I-2mu)>aJpha>0  results  in  "leaky  LMS") 

•  OUTPUTS: 

•  y  -  lxT  vector  of  LMS  filter  outputs 

•  err  -  lxT  vector  of  LMS  error  outputs 

•  wts  -  MxN  vector  of  the  final  converged  filter  weights 

•  ALGORITHM: 

•  The  Leaky/Std  LMS  algorithm  is  defined  as  follows  at  time  k, 

•  1.  y(k)  *  sum  (m*  1 M  sum  (i»  1)  **  N  conj(wts(mj))*x(mJt-i) 

•  2.  err(k)  *  d(k)  -  y(k) 

•  3.Ex(k)-  W‘Ex(k-I)  +  (l-W)-sum  (m*1>  |x(k)|  "2 

•  4.  mu'  «  (2“mu)  /  (eo*  +  N*Ex(k)) 

•  5.  wts(k+ 1  Im^))  *  (1- mu’ *  alpha  )*wts(k|m,p)  + 

•  mu’4conj(erT(k))*x(mJt-p) 

•  For  alpha  «0,  the  Leaky  alogirthm  reduces  to  the  Std  LMS 

•  algorithm.  The  filter  is  configured  as  an  Adaptive  Noise 

•  Canceller. 

•  REFERENCE: 

•  S.  Haykin,  "Adaptive  Filter  Theory,",  Prentice  Hall.  1991. 

•  written  by:  Richard  North  MJ2 

. . . . . . . . * . * . ) 

#mdude  <math.h> 

♦include  <stdio.h> 

♦include  "cmexh" 


%  function  [Py.y]  ■  td_cbf(  X,  az.  el,  P,  dt,  c) 


♦define  MAX  N  151  /*  Maximum  filter  Length  for  Momentum  LMS  V 
♦define  EPS  'D.OOOOOI  r  Minimum  number  ■  1x10-6  •/ 


mmmmmm  mmmmmmmmmmmmmmmm  mmmmmm  f  input  arguments  V 

%  initialize  variables  ♦define  X_1N  prhs[0] 

%mmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm  ♦define  D  IN  prhslll 

fM,T]  -  size(X);  ♦define  NJN  prhsf2l 

B  *  1/dt;  %  sampling  frequency  ♦define  mu  IN  prhspl 

n_a2  -  length(az);  ♦define  alpha JN  prhs]4] 

n”el  ■  length(el); 

NIooks  -  n_az*n  el;  r  output  arguments  V 

%W  ■  ones(M,l  f;  %  unity  sensor  weighting  ♦define  Ypred_OUT  plhs{0] 


Figure  A 2.  Continued. 
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#dcfine  E n  OUT  plhslll 
#definc  Wtf  OUT  plhs[2) 


user  fcn(nlhs,plhv>rhs,prhs) 
int  rilhs,nrhs; 

Matrix  *ptiu[),  *prhs|]; 


/•  Matlab  variables  7 
unsigned  int  T; 
double  N,mu,alpba; 

double  •X_re,*D_re,*Ypred  re,*ErT_re,*Wts  re; 
double  ,X~im,,D”  im/Ypretl  im.’Err  lm.'Wts  im; 
double  Wts_temp[MAX_N].Wts_old|  MAX_N);  " 


r  mi  sc  local  variables  7 

int  t,m,n,dl42>M; 

double  Ex,  normmu.W* 0.999; 


/•  check  number  of  arguments  7 

if  (nrhs  !■=  5)  mex  error(”'lms’  requires  5  inputs.”); 

if  ((nlhs  >  2)  &&  Tnlhs  !=  3))  mex^error<”’lms’  requires  3  outputs.”); 

if  (nlhs  <  2)  mex_error(’”lms‘  requires  at  least  2  outputs.”); 


/•  set  variables  and  allocate  memory  7 
N**NIN->pr;  /*#  of  filter  taps  7 

if  (N  >  MAX  N)  mex_error(”’lms'  N  >  MAX  N  need  to  recompile  lms.c”); 
mu  *  *mu  1ST-  >pr;  ~  f*  value  of  step  sixe  *7 
alpha  *=  *aTpha_IN->pr;  /*  value  of  leak  parameter  7 
T  ■  X  IN  -  >nT  I*  #  of  columns  equals  T  7 
M  «  >m;  /•  #  of  rows  equals  M  (no.  of  channels)  7 

X  re  -  X  lN->pr; 

X_im  =  X~  IN — >pi; 

D~re  •=  DTN->pr; 

D~un  “  D^_IN- >pi, 


if  ((X  un  -  -  NULL)  &&  (D_im  !  *  NULL))  /•  x  is  real  &  d  is  complex  7 
mex_error(”x  is  real  &  d  is  complex"); 
if  ((X~im  !«  NULL)  &&  (D_im  ■  *  NULL))  r  x  is  complex  &  d  is  real  7 
mex~error(”x  is  complex  &~d  is  real”); 

if  |X_im  --  NULL)  }'  real  7 

Ypred  OUT  «  create_matrix(l.T,REAL); 

Err  OTJT  «  create  mitnx(l,T,REAL); 

Wts^OUT  ■  CTeate_matrix(M,(int)N,REAL); 

Ypred_re  *  Ypred  OUT->pr; 

Err  re- «  Err  OUT-  >pr; 

Wts"  re  -  Wts  OUT-  >pr; 

else  r  complex  7 

Ypred  OUT  -  create  mat  nx(l,T.  COMPLEX); 

Err  OtJT  *  create  matrix(l,T,COMPLEX); 

Wts' OUT  =  create' matrix(M,(int)N.COMPLEX); 

Ypred  re  ®  Ypred  OUT->pr; 

Ypred fim  ■  Ypred-  OUT->pi; 

Err_re  ®  En_OUT-  >pr, 

Err'im  *  Erf  OUT->pi; 

Wts'  re  -  Wts'  OUT->pr; 

Wtsjm  *  Wts“OUT->pi; 


if  jX_im  --  NULL)  r  real  7 

/•  NOTE:  Matlab  stores  by  columns  so  that  X(ij)*X[i+Nrow$*j)  7 
r  initialize  real  variables  7 
for  (<®0;t<T;  t+  +  ) 

Ypred  re[t)  *  0.0; 
for  (m*0;  m<M;  m++) 
for  (n*0;  n<N;  n+  +  ) 

Wts  re|m+M*n)  ■  0.0; 

Ex -Of 

for  (t»0;  t<(int)N-l;  t+  +  ) 

for(m*0;m<M;  m++) 

Ex  -  W'Ex  +  (l-W)'(X_re[m+M,tj,X_re[m+M*l]); 


/*  filter  real  data  with  real  LMS  7 
for  (t«(int)N-l;  t<T;  t++) 

for  (m«0;  m<M;  m++) 

for  (n®0;  n<(int)N;  n++) 

Ypred_re[t]  +-  Wts_re(m+M*n)*X_refm+ M*(t-n)J; 

Ex  -  W'Ex  +  (I -W)*(X_re(m+M*t],X_re(m+M*t}); 

tf7&*>tps)-r‘,')"Ypreil-rc,‘1; 

norm  mu  -  (2*mu)  /  (N*Ex); 
else 

normjnu  -  (2*mu)  /  (N*EPS); 

for  (m»0;  m<M;  m++) 

{ 

for  (n*0;  n<(int)N;  n+  +  ) 
r  Leaky/Std  LMS  7 

Wts_re|m+M*n]  ■  (1.0  -  norm_mu*alpha)*Wts_re(m+M*nj  -f 


norm^mu  “  En_re  [t )  *  X_re  J  m + M  *  1 1  -  n )): 

) 

} 

) 

else  /•  complex  7 

r  NOTE:  Matlab  stores  by  columns  so  that  X(i,j)=X(i+Nrows"j]  7 
r  initialize  complex  variables  7 
for  (t=0;  i<T;  t++) 

( 

Ypred_re(t]  *  0.0; 

Yprcd_un[t)  *  0.0; 

for  (m*0;  m<M;  m+  +  ) 
for  (n*0;n<N;n+  +  ) 

< 

Wts  re[m+M*n)  *=  0.0; 

Wts~im(m+M*n]  «  0.0; 

Ex  *  0; 

for  (t*0;  t<(int)N-l;  t+  +  ) 

for  (m®0;  m<M;  m++) 

dl  -  m+MT; 

Ex  -  W’Ex  +  (1-W)*(X  re(dl)’X  re[dl)  +  X  im(dl)*X  im(dl)); 


r  filter  complex  data  with  complex  LMS  7 
for  (t*(int)N-l;  t<T;  t+  +  ) 

for  (m®0;  m<M;  m++) 

< 

for  (n*0;  n<(int)N;  n++) 

dl  *  m+M*n; 
d2  =  m+M"(t-n); 

Ypred  reft)  +«=  Wts  re|dl)*X  re(d2]  +  Wts  im[dl)-X  imld2); 
Ypred'im[t)  +*  Wts_re(dl)“>T^imIdi]  -  Wts_im(dI]*X_re[d2); 

dl  *=  m+M*t; 

Ex  *  W*Ex  +  (l-W)-(X_re[dl)'X_re(dl]  +  X_im[dl)*X^im[dl]); 

ErT  reft]  «  D  reftj- Ypred  reft); 

ErTim[t^=^Iir  imft)  -  Y  pre  J_im(t]; 

norm  mu  ■  2*mu/(N*Ex); 
else 

normjnu  *  2*mu  /  (N’EPS); 

for  (m*0;  m<M;  m+  +  ) 

for  (n®0;  n<(int)N;  n++) 

dl  -  m+M'n; 
d2  «  m+M#(t-n); 
r  Leaky/Std  LMS  7 

Wts  refdlj  ®  (1.0  -  norm  mu*atpha)*Wts  rejdll  + 

'norm  mu*(EjT_re[t)’X_re[d2j  +  Ett  un[tl*X  im[d2)); 
Wtsjmfdf]  «  (1.0- norm  mu^alphaJ^Wts  tm[dir+ 

'Dorm_mu*(ErT_re[t]*X_im[dz)  -  ErT_un[t]*X_re(d2]); 

) 

} 


gsc.m 


function  [Py.Y]  ■  gsc<  X.  az,  el.  P.  dt.  c.  N.  mu,  B) 

%  GSC  -  Generalized  Sidelobe  Canceller  (Linear  Constrained  Adaptive 
%  Beamfonner) 

% 

%  INPUT: 

%  X  -  MxT  array  of  T  time  samples  of  the  M  sensor  outputs 
%  az  -  beamforming  azimuth  angles  (degrees) 

%  el  -  beamforming  elevation  angles  (degrees) 

%  P  -  3xM  a  nay  of  the  M  sensor  positions  (ft)  (m) 

%  dt  -  sample  increment  (sec) 

%  c  -  speed  of  sound  in  propagation  media  (ft/sec)  (m/sec) 

%  N  -  channel  filter  length  (same  for  all  elements) 

%  mu  -  LMS  step  size 
%  B  -  ixM  blocking  matrix 

%  (where  J  <  M  is  the  number  of  sidelobe  cancellers) 

% 

%  OUTPUT: 

%  Py  -  array  of  beamfonner  power  output;  Py(kJ)  is  the  output 
%  at  azimuth  angle  az(l)  and  elevation  angle  el(k). 

%  Y  -  NlooksxT  anay  of  T  time  samples  of  the  beamfonner  output 
% 

%  ALGORITHM  . 

%  As  discussed  in  (1)  A  |2],  the  Linear  Constrained  Adaptive  Beamfonner 
%  is  implemented  iteratively  as  a  Generalized  Sidelobe  Canceller  (GSC). 
%  It  consists  of  a  nonadaptive  time  domain  classical  beamformer  followed 
%  by  a  FIR  filter.  To  cancel  interferers  which  bleed  in  through  the  CBF 
%  side  lobes,  a  sidelobe  canceller  (consulting  of  a  blocking  matrix  and 
%  a  multichannel  adaptive  filter)  adaptively  filters  the  interferes  in 
%  space  and  time  to  subtract  from  the  CBF  output. 

% 

%  NOTE; 

%  Two  possibilities  for  B  are  M 1 

%  1.  col  •  I  1  zeros(l.M-z));  2.  col  ■  (  1  zen»(l.M-3)); 
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%  row  *  ( 1  -1  zeros(l,M-2)];  row  *  [  1  -2  1  zeros(l,M-3)]; 

%  B  «  loepliu(coUow); 

% 

%  REFERENCES: 

%  [1]  L.  Griffiths,  C.W.  Jim, “An  Alternative  Approach  to  Linear  Constrained 
%  Adaptive  Beamfomung" 

%  IEEE  Trans,  Antennas  and  Prop.,  vol.  AP-30.  pp.  27  -  34. 1982. 

%  [2j  R.  North.  "Broadband  Conventional  Beamfomung  Incorporating 
%  Adaptive  Equalization NRaD  Tech.  Note  1992 

% 

%  function  [Py.Y]  «  gsc(  X.  az,  el,  P.  dt,  c,  N,  mu.  B) 

*.«««■»= BBS*.*.... 


%  initialize  variables 

[M.T]  *  size(X); 

J.M11  -  size(B); 

If  ((Ml  M)  |  (J  >«  M)) 
error('Blocking  Matrix.  B,  wrong  size!'); 
end; 

fs  *  1/dt;  %  sampling  frequency 
n_az  *  length(az); 
n“el  «=  length(el); 

Nlooks  *  n_az*n  el; 

%W  *  ones(M,l7;  %  imity  sensor  weighting  for  CBF  (see  also  chebwin) 
Y  «  zeros(N1ooks,T); 

Yc  -  zeros(l,T); 

Ya  *  zeros(l,T); 

Py  *  zeros(l, Nlooks), 

Xa  «  zeros(M.T); 

Xp  ®  zcros(J.T); 

%  Generate  steering  time  delays  for  all  directions  of  look 


A  *  steer(  az ,  el,  P ,  c ,  1 ,  Time');  %  delays  in  seconds  (freq  not  used) 
TAU  *  round(A'fs);  %  delays  in  iterations 


%  align  data  matrix  and  beamform  in  each  look  direction 

for  n=l:NIooks 
%  align  all  element  data 
Xa  «  xalign(X,TAU(.\n)); 

%  form  sidelobe  canceller  input  matrix 
Xp  -  B  *  Xa; 


%  weight  all  element  data  for  TD-CBF 
%Xa  «  diag(W)*Xa; 

%  form  TD-CBF  output 
Yc  *  sum(Xa)/M; 

%  filter  TD-CBF  output 
%Yc  ■  fiIter(b,a,Yc); 

%  multichannel  adaptive  filter  Xa  matrix 
( Ya,  Y(n,:  ),W]  *  mch_nlms(Xp,Yc,N  jnu.O); 
end; 

Py  -  mean(Y\*conj(Y’)); 

%  tidy  up  and  reshape  the  output; 


%  we  know  Py  must  be  real  ... 

Py  -  real(Py); 

%  ...  and  return  v  with  the  rows  containing  the  azimuth  for  \  fixed 
%  elevation  angle. 

Py  ■  reshape(  Py ,  n_az .  n_el )'; 

inch  af.m 


function  fPy.Y]  »  mch_af(  X  az,  cl,  P,  dt,  c,  N,  mu,  d) 


%  MCH  AF  -  Multichannel  Adaptive  Filter  with  presteering  (data  alignment) 
% 

%  INPUT: 

%  X  -  MxT  array  of  T  time  samples  of  the  M  sensor  outputs 
%  az  -  beamfomung  azimuth  angles  (degrees) 

%  el  -  beamforming  elevation  angles  (degrees! 

%  P  -  3xM  troy  of  the  M  sensor  positions  (ft)  (m) 

%  dt  -  sample  increment  (sec) 

%  c  -  speed  of  sound  in  propagation  media  (ft/secl  (m/sec) 

%  N  -  channel  filter  length  (same  for  all  elements) 

%  mu  -  LMS  step  size 
%  d  -  IxT  vector  of  desired  values 
% 

%  OUTPUT: 

%  Py  -  array  of  beamfbrmer  power  output:  Py(k,l!  is  the  output 
%  at  azimuth  angle  az(l)  and  elevauon  angle  el(k). 

%  Y  -  NlookixT  array  of  T  time  samples  of  the  beamformer  output 
% 

%  REFERENCES: 

%  [1]  R.  North,  "Broadband  Conventional  Beamforming  Incorporating 


%  Adaptive  Equalization,"  NRaD  Tech.  Note  1992 

% 

%  function  (Py.Y)  »  mch_af{  X  az.  el.  P,  dt.  c,  N,  mu.  d) 


%  initialize  variables 
(M,TJ  *  sue(X). 

fs  *  1/dt;  %  sampling  frequency 
n_az  «  iength(az); 
n”el  *  length(el); 

NIooks  =  n_az*n  el; 

%W  =  oncs(M,lf;  %  unity  sensor  weighting  for  CBF  (see  also  chebwin) 
Py  =  zeros(l, Nlooks); 

Xa  *  zeros(M.T); 

%  Generate  steering  tune  delays  for  all  directions  of  look 

A  *  steer(  az  .  el,  P ,  c .  1  ,  Time*);  %  delays  in  seconds  (freq  not  used) 

TAU  s  round(A'fs);  %  delays  in  iterations 


% 


Shift  Delays  so  that  the  Filtering  is  centered  giving  equal  +/-  looks 


TAU  *  TAU  -  N/2*ones(TAU); 
MinTAU  =  min(min(TAU)l; 

MaxTAU  *  max(max(TAU)); 
if  MinTAU  <  0 
Tmax  *  T+ MinTAU; 
else 

Tmax  *  T; 
end 

if  MaxTAU  >0 
Tmin  *  MaxTAU+1; 
else 

Tmin  *  1; 
end 

Y  .  zeros(Nlooks3ength(Tmin:Tmax)); 


%  shift  all  delays  by  N/2 


%mmmmmmmmmm ....... b . b . b . b b . b . - - ..... 

%  align  data  matrix  and  beamform  in  each  look  direction 

for  n=*l:Nk>oks 
%  align  all  element  data 
Xa  «  xalign(XTAU(:,n)); 

%  weight  all  element  data 
%Xa  -  diag(W)*Xa, 

%  multichannel  adaptive  filler  nonzero  Xa  data 
(Y(n,:XE,W).mch_nlms(Xa(:,Tntin:Tmax)^(Tmin:Tmax),Njnu.O); 
end; 


Py  ■=  mcan(Y’."conj(Y’)); 

%  tidy  up  and  reshape  the  output; 


%  we  know  Py  must  be  real ... 

Py  *  real(Py); 

%  ...  and  return  y  with  the  rows  containing  the  azimuth  for  a  fixed 
%  elevation  angle. 

Py  .  reshape(  ry ,  n_az ,  n_el )’; 


%. 


arraypat.m 


function  W  .  arraypat(  az.el.P.c.fO.u.w) 


%  ARRAYPAT  -  computes  the  amy  pattern  assuming  unity  shading 
%  and  plane-wave  inputs. 

% 

%  INPUT: 

%  az  -  vector  of  azimuth  angles  to  sweep  (deg) 

%  el  -  vector  of  elevation  angles  to  sweep  (deg) 

%  P  -  3xM  array  of  the  M  sensor  positions 
%  c  -  speed  of  signal  propagation 
%  10  -  frequency  to  beamform 
%  u  -  array  steered  direction  vector  (az  el]  (deg) 

%  w  -  lxM  vector  of  element  shading  coefficients 
%  (if  w  is  not  entered,  then  unity  shading  is  assumed) 

% 

%  OUTPUT;  let  M  denote  the  number  of  sensors;  set 
%  n_az  *  length(az) 

%  n“el  *  length(el) 

%  Then  the  output  army  pattern  W  is  n_azxn  el  where  each  column 
%  has  the  for 

%  T 

%  W(:Jt)  -  sum  (m«l)  ~  (m*M>  w_m*[exp(  i  2  pi  fOP  (A(^k)-u_sVc)] 
% 

%  where  u_s  (3x1)  is  the xyz  coordinates  of  the  array  steered 
%  direction, 

% 

%  U_1  -  |  coKu(l))coi(u(2))  | 

%  |  «m(u(I))CM(u(2))  | 

*  I  | 

% 

9b  and  A(:Jt)  (3xn_az*n_el)is  the  xyz  matrix  of  assumed  incident 
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%  narrowband  complex  exponential  sources 
% 

%  A(l,k)  -  |  cos(az(l))  cos(el(k))  | 

%  I  sm(a2(l))cos(el(k))  ] 

%  |  sm(el(k))  | 

% 

%  written:  Rich  North  5-19-92 
%  W«arTaypat(az,el.P.c.f0.u)  or 
%  W«anaypat(az,el.P,c,f0.u.w) 


n_az  *  length(az); 
n^cl  «  length!  el); 
[nrow.Mj  «  size(P); 
if  nargrn  <  7 
w  ■  ones(l,M); 
end 


%  use  steer.m  to  generate  time  delays 
A  *  steer(  az ,  el .  P  .  c  .  fO ,  Time’ ); 
u_s  =  steer(  u(l) ,  u(2) ,  P  ,  c  ,  fO  ,  Time’ ); 

T  *  exp(i  *  2  *  pi  *  fD  *  (A-u_s*ones(l,n_az*n_el))); 

%  sum  all  shaded  elements  contribution  to  the  spatial  component  of  total 
%  array  output  (time -domain) 

W  *  onesfl,n_a2*n  el); 

[nrow^icol]  »"size(w); 
if  nrow  m  m  l  %  w  ■  row  vector 
W  -  w*T/M; 

elscif  ncol  «  *  1  %  w  -  column  vector 
W  *  w"T/M; 
else 

crror(’...Error  w  should  be  a  row  vector  lxM’); 
end 

W  -  reshape(  W ,  n_az ,  n_el )’; 

bfclass.m  (modified  slightly  from  [6]) 


function  y  ■  bfclass(  X,  R,  fl),  az,  el.  P.  dt,  c) 

—  ■  —  ■  —  —  ■ 

%  BFCLASS  -  beamforming  using  the  classical  estimation  procedure 

% 

%  INPUT: 

%  X  -  MxN  array  of  N  time  samples  of  the  M  sensor  outputs 
%  R  -  Estimate  the  MxM  cross -spectral  density  matrix  at  fO: 

%  fD  -  frequency  to  operate  beamformer  (Hz) 

%  n_sg  -  number  of  points  per  block  average  segment  -  >  must  divide  N 
%  az  -  beamforming  azimuth  angles  (degrees) 

%  el  -  beamformins  elevation  angles  (degrees! 

%  P  -  3xM  array  of  the  M  sensor  positions  (ft) 

%  dt  -  sample  increment  (sec) 

%  c  —  speed  of  sound  in  propagation  media  (ft/sec) 

% 

%  OUTPUT: 

%  y  -  array  of  classical  beamformer  output;  y(kj)  is  the  output 
%  at  azimuth  angle  az(l)  and  elevation  angle  el(k). 

% 

%  ALGORITHM:  As  discussed  in  (1).  the  classical  beamformer  has  output 
%  H 
%  y*a  Ra 

% 

%  for  steering  vector 
%  H 

%  a  «  [exp(i  2pifDp_(m)  u/c)] 

% 

%  where  p_fm}  is  the  position  vector  of  the  mth  sensor,  u  is  the  look 
%  direction. 

% 

%  REMARKS: 

%  (1 ) |  to  be  am  form  at  a  desired  frequency  fO  (Hz)  requires  that  the  sample 
%  increment  dt  (sec)  first  satisfy 
%  I0<  l/2dt 

%  (2)  y’  is  the  conjugate  transpose  of  the  vector  y 

% 

%  REFERENCES: 

%  (11  Don  H.  Johnson. 

%  'The  Application  of  Spectra)  Estimation  Methods  to  Bearing 
%  Estimation  Problems  , 

%  Proceedings  of  the  IEEE.  Volume  70,  Number  9,  pages  1018- 1028, 1982 
% 

%  function  y  -  bfclass(  X.  R,  fO,  az.  el,  P,  dt,  c) 


%  Get  the  array  A  of  steering  vectors 
n_az  ■  length(az); 
n~el  ■  length!  el); 

A“  -  steer(  az .  el,  P ,  c .  fO ); 


%  Be  am  form  using  the  correlation  matrix  R  at  the  specified  azimuth  angles 
%  and  the  assumed  input  frequency  fD. 
y  -  zeros(  1  ,n  az*n  el); 
fork*l:(n  az*n  ef) 
y(k)-AfrtrrR’A(:,k); 
end 


%  tidy  up  and  reshape  the  output; 
%  we  know  y  must  be  real  . 
y  -  real(y); 


%  ...  and  return  v  with  the  rows  containing  the  azimuth  tor  a  fixed 
%  elevation  angle, 
y  *  reshape!  y  ,  n_az  .  n_el )'; 

mvdr.m  (modified  slightly  from  [6]) 


function  (y.w)  *  mvdr(  X.  R.  fO,  az.  el,  P.  dt,  c) 

%«*****  =  «-*  =  **■■*  ■  * 

%  MVDR  -  Minimum- Variance  Distortionless  Response  beamforming 

% 

%  INPUT; 

%  X  -  MxN  array  of  N  time  samples  of  the  M  sensor  outputs 
%  R  -  Estimate  of  the  MxM  cross -spectral  density  matrix  at  fO: 

fl)  -  frequency  to  operate  beamformer  (Hz) 

%  n_sg  -  number  of  points  in  each  block  averaging  segment  -  >  must  divide  N 
%  az  -  beamforming  azimuth  angles  (degrees) 

%  el  -  beamforming  elevation  angles  (degrees) 

%  P  -  3xM  array  of  the  M  sensor  positions  (ft) 

%  dt  —  sample  increment  (sec) 

%  c  -  speed  of  sound  in  propagation  media  (ft/sec) 

%  OUTPUT: 

%  y  -  MVDR  beamformer  output;  y(kj)  is  the  output  at  azimuth  angle 
%  az(l)  and  elevation  angle  ei(k). 

%  w  -  (n_az*n_ellxM  array  of  row  vectors  containing  the 
%  element  weignts  for  each  azimuth  and  elevation  angle 
%  (only  outputed  if  nargout  *  2) 

% 

%  ALGORITHM:  As  analyzed  in  [2],  the  MVDR  beamformer  has  output 
%  H  H  + 

%  y  -  w  R  w  -  l/(a  R  a)  (1) 

%  w  «*  (R  a)/(a  R  a)  (2) 

% 

%  for  steering  vector 
%  H 

%  a  *  (exp(i  2  pi  f  p_{m)  u/c  )]  (m*1.2, ....  M) 

% 

%  + 

%  R 

%  is  the  Moore -Penrose  inverse  R  PROVIDED  the  steering  vector  is 
%  orthogonal  to  the  null  space  of  R.  Otherwise,  if  the  steering 
%  vector  has  a  non -zero  component  in  the  null  space  of  R,  then 
%  the  MVDR  output  is  y  *  0. 

% 

%  WARNING: 

%  This  function  MVDR  assumes  that  the  estimated  CSD  matrix  R  is  full 
%  rank  and  therefore  will  calculate  y  as  in  equation  ( 1 ). 

% 

%  REMARKS: 

%  (1)  to  beamform  at  a  desired  frequency  fO  (Hz)  requires  that  the  sample 
%  increment  dt  (sec)  first  satisfy 
%  fO  <  l/2dt 

%  f  2)  y’  is  the  conjugate  transpose  of  the  vector  y 
%  (3)  to  obtain  a  hill  rank  CSD  estimate,  it  is  necessary  to  use  at  least 
%  M  block  averages  in  function  "sfecsd". 

% 

%  REFERENCES: 

%  (11  Don  H.  Johnson, 

%  'The  Application  of  Spectral  Estimation  Methods  to  Bearing 
%  Estimation  Problems’’, 

%  Proceedings  of  the  IEEE,  Volume  70.  Number  9,  pages  1018- 1028. 1982 
%  (2)  Jeffrey  Speiser 

%  “MVDR  Fundamentals  and  Computation” 

%  preprint.  Code  635,  NOSC,  San  Diego.  CA,  92152-  5000 
% 

%  y  *  mvdr(  X,  R.  fO.  az,  el,  P.  dt,  c)  or 
%  (y.w)  -  mvdr(  X.  R,  fO.  az.  el.  P,  dt,  c) 


%  Calculate  the  inverse  of  the  CSD  matrix:  it  is  assumed  that  R  is  full 
%  rank  and  therefore  the  Moore -Penrose  inverse  and  inv(R)  are  the  same. 
R_pinv  ■  pinv(R); 

%  Get  the  array  of  steering  vectors. 

%  On  output,  the  matrix  A  will  be  M  x  (  n_az  n_el ). 
n_az  *  length(az); 
tfel  ■  length(el); 

/T*  steeifaz ,  el ,  P ,  c .  (0  ); 


%  Beamform  using  the  correlation  matrix  R  at  the  specified  azimuth  angles 
%  and  the  assumed  input  frequency  fO 
y  -  zeros(  1,  (n_az*n_el) ); 
foTk«l:(n  az*n  el) 
y(k)  -  AC-XY  '■Rjxot  •  A(U); 
end 


%  we  know  y  must  be  real ... 
y  -  I  yrealfy); 


if  nargout  -«  2 

%  generate  the  weights  for  each  azimuth  and  elevation  steering 
|nrow,M]  -  sizefP); 
w  ■  zero*((n  az*n  el),M); 
w  ■  (R_pmvn  A).r.*  (y**oncs(l,M)); 
end 


%  tidy  up  and  reshape  y: 

%  ...  the  k-  th  row  of  y  contains  the  beamformer  output  at  elevation 
%  ei(k)  and  all  the  aziumth  angles, 
y  -  reshape!  y ,  n_az ,  n^el )’; 
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Appendix  B 

TIME-DOMAIN  CBF  OF  A  LINEAR  ARRAY  WITH  UNIFORMLY 

SPACED  ELEMENTS 


As  an  example  of  time-domain  CBF  (and  a  more  detailed  CBF  review),  consider  a 
linear  array  with  uniform  element  spacing,  d,  and  a  single  source  such  that 

x m  =  (m  d)  z  m  =  0, 1,2,..., A/  -  1 

?!  =  (cos (f>1  cos0j)  x  +  (cos0!  sin0j)  y  +  (sin0j)  z 

u  =  (cos0  cos0)  x  +  (cos <p  sin 0)  y  +  (sin0j  z 

It  is  easy  to  observe  that  the  delay  required  for  coherent  summing  across  the  array  requires 
u  =  Sj,  or  a  delay  equal  to 


T„(U)  = 


=  round 


c  c  (sec) 

m  d  sin <p,  ,  1  ..  . 

- c - 1  f sample  I  (iterations) 


(B.l) 

(B.2) 


Figure  B.l  shows  the  uniformly  spaced,  linear,  vertical  array  (ULVA)  and  a  single  source. 
When  the  array  is  not  steered  directly  at  the  source,  Eq.  NO  TAG  becomes 


m  =  0 
M-  1 


])4 

n  =  0  V  L 


y(n,u)  =  £  wm  njn  -  round  —  ^  fsample 
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Figure  B.l.  Single  source  ULVA. 

To  proceed  analytically,  we  must  make  some  assumptions  about  the  source  model. 
Assume  a  sampled,  narrowband  complex  exponential  source,  s^n)  —  A j  where 

f\  ~  f\  /  fsample  *s  t*ie  normalized  source  frequency.  Then  Eq.  (B.3)  can  be  rewritten  as 
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y(n,u)  =  X  wm  njn  -  round  ^  sintp  fsampleU  + 

m  =  0  '  ' 
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And  the  power  spectrum  of  Eq.  (B.4)  for  stationary  white  noise,  nm(n),  is  [6] 
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Note  that  the  second  term  of  Eq.  (B.4)  is  composed  of  a  time-varying  component  (left 
term)  and  a  spatially-varying,  time-invariant  component  (right  term).  The  spatially-varying 

f 

term  is  a  function  of  wavelength,  X  =  element  spacing,  and  the  direction  of  the  source 
relative  to  the  array.  The  spatially-varying  component  is  usually  written  in  terms  of  the 
wavenumber,  K  =  Thus,  we  find  the  continuous-time  spatial  component  of  the  array  is 

A 

equal  to 

M—  1 

W(<t>)  =  X  wm  e-jK>md(sin<f>-sin<p,)  (B.6) 


m  =0 

and  the  discrete-time  spatial  component  of  the  array  equal  to 
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Eq.  (B.6)  can  be  viewed  as  the  spatial  transfer  function  of  the  array,  since  we  can  let 
z  as  e;K,</(s»n0-sm0,)  an(j  rewrite  Eq.  (B.6)  as  a  (right-sided  i.e.,  appropriate  for  causal 

M-l 

sequences)  Z  transform,  W(<f>)  =  X  Wm  z~m-  Equation  (B.7)  can  be  thought  of  as  the 

m  =  0 

sampled  spatial  transfer  function.  Thus,  the  number  of  zeros  defining  the  spatial  component  of 
the  array  is  equal  to  the  number  of  elements,  M.  It  is  apparent  that  the  shading  coefficients 
behave  similarly  to  standard  FFT  weighting  [2]  for  the  uniformly  spaced  linear  array. 
However,  because  of  the  nonlinear  relationship  between  W(<p)  and  <pv  the  3-dB  beamwidth  of 
the  main  lobe  is  minimum  at  <f>  =  0  degrees  (broadside)  and  maximum  at  <p  =  ±  90  degrees 
(end-fire). 


The  spatial  filtering  characteristics  (called  the  array  pattern  for  the  general  array 
configuration  case)  of  an  8 -element  ULVA  with  half-wavelength  spacing  and  no  shading  are 
demonstrated  in  Figures  B.2  and  B.3  for  the  array  steered  in  the  broadside  direction  and  the 
near  end-fire  direction.  It  is  easily  seen  from  the  figures  that  the  main  lobe’s  spatial  width 
increases  as  the  steering  angle  approaches  +/—  90  degrees.  It  was  found  in  [1]  that  the  main 
lobe  beamwidth  is  approximately  equal  to 

BW  «  49.6  ^  sec{<p)  (degrees)  (B.8) 
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for  steering  angle  0.  The  effects  of  sampling  the  continuous  data  are  seen  in  figures  B.4  and 
B.5  for  fsanpie  =  250  Hz  and  fsamp[e  =  1500  Hz,  respectively.  The  delineation  of  the  array 
pattern  is  seen  to  improve  as  the  sampling  frequency  increases.  Finally,  figures  B.6  and  B.7 
show  the  ULVA  pattern  for  detecting  frequencies  for  which  the  array  is  not  exactly  tuned,  i.e., 

j  *  0.5.  Referring  back  to  the  Z-transform  description  of  the  array  pattern,  it  is  seen  that  the 

Z-plane  unit  circle  is  traversed  exactly  2  y  times.  Thus,  for  y  <  0.5,  the  source  frequency  is 
less  than  the  designed  frequency  and  the  array  pattern  only  traverses  part  of  the  Z-plane  unit 
circle.  For^  >  0.5,  the  source  frequency  is  larger  than  the  designed  frequency  and  the  array 

pattern  traverses  the  Z-plane  unit  circle  more  than  once,  resulting  in  a  type  of  aliasing.  This 
generates  spurious  lobes  called  “grating  lobes”  which  create  an  ambiguity  as  to  the  true 

direction  of  arrival.  Only  for  j  =  0.5  is  the  Z-plane  unit  circle  traversed  exactly  once. 

The  time-domain  CBF  can  be  demonstrated  for  two  narrowband  inputs  separated  by 
only  5  degrees.  A  40-element  ULVA  is  used  with  low  SNRs  for  both  inputs.  Figure  B.8  plots 
the  beampattem  averaged  over  100  iterations  for  <j)l  =  0  degrees,  =  25  Hz,  SNR}  =  0  dB 
and  <p2  =  -  5  degrees,  f2  -  45  Hz,  SNR2  =  -  5  dB.  From  figure  B.8  one  can  see  the 
excellent  results  which  can  be  obtained  with  time-domain  CBF.  Appendix  A  lists  the  MATLAB 
code  used  to  generate  the  ULVA  data. 
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Figure  B.2.  ULVA  pattern  for  j  =  0.5,  <pm  =  0  degrees,  N  =  8,  fsompie  =  00 . 
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Figure  B.4.  ULVA  pattern  for  j  =  0.5,  0  =  0  degrees,  M  =  8,  /Mwpfe  =  250  Hz. 
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sin(phil)  linear  vertical  array 


Figure  B.5.  ULYA  pattern  for  ~  =  0.5,  (j>  =  0  degrees,  M  =  8,  fsampie  =  1500  Hz. 


-1  -0.8  -0.6  -0.4  -0.2  0  02  0.4  0.6  0.8  1 

«in{phil)  linear  vertical  imy 


Figure  B.6.  ULYA  pattern  for  j  =  0.25,  <p  =  45  degrees,  M  =  8,  /sampfe  =  ® . 
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